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Chapter  1 
Introduction 


An  artificial  medium  is  basically  a  macroscopic  analog  of  a  real  medium,  and  typically 
consists  of  a  large  number  of  scattering  objects  distributed  (more  or  less)  uniformly 
in  some  host  or  background  medium.  The  scattering  objects  affect  the  behavior  of 
electric  and  magnetic  fields  inside  the  artificial  medium.  For  example,  a  general  plane 
wave  in  an  artificial  medium  propagates  with  an  effective  wavenumber  different  than 
that  of  the  host  material  or  the  scattering  objects.  As  a  result,  artificial  media  can 
be  characterized  by  an  effective  permittivity  and  effective  permeability.  In  general, 
artificial  media  are  anisotropic,  and  the  constitutive  parameters  are  tensor  quantities. 

When  a  plane  wave  propagates  through  an  artificial  medium,  currents  are  induced 
in  (or  on)  the  scattering  objects.  As  shown  in  Figure  1.1,  these  currents  can  be 
viewed  as  macroscopic  current  moments,  analogous  to  the  microscopic  dipole  moments 
induced  in  the  molecules  of  an  actual  dielectric  [1].  The  effect  of  the  macroscopic 
current  moments  is  to  produce  a  net  electric  and  magnetic  current  moment  per  unit 
volume,  and  thus  the  artificial  medium  has  some  complex  effective  permittivity  and 
permeability  different  from  the  host  medium. 

The  complex  effective  constitutive  parameters  of  the  artificial  medium  are  a  func¬ 
tion  of  frequency,  the  electrical  size,  shape,  spacing,  and  orientation  of  the  scattering 
objects,  and  the  constitutive  parameters  of  both  the  host  medium  and  the  scattering 
objects.  Also,  in  contrast  to  real  media,  the  constitutive  parameters  can  be  a  function 
of  the  direction  of  propagation  through  the  artificial  medium.  By  properly  choosing 
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Artificial  Dielectric 


(b) 


Figure  1.1:  The  artificial  dielectric  model. 
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the  geometry  and  composition,  it  may  be  possible  to  design  an  artificial  medium  of 
desired  permittivity,  permeability  and  loss  tangent. 

This  dissertation  presents  an  integral  equation  and  periodic  method  of  moments 
(PMM)  [2,  3,  4]  solution  to  the  problem  of  determining  the  effective  permittivity 
and  permeability  of  an  artificial  medium.  The  artificial  medium  is  composed  of  a 
3D  periodic  array  of  identical  arbitrarily-shaped  thin  conductive  or  dielectric  wire 
objects. 

The  solution  proceeds  as  follows.  First,  an  integral  equation  is  formulated  for  a 
plane  wave  of  unknown  wavenumber  propagating  in  an  artificial  medium  of  infinite 
extent  in  all  three  dimensions.  Next,  this  integral  equation  is  solved  by  the  PMM, 
yielding  the  complex  effective  wavenumber  of  the  plane  wave,  the  eigenfunction  cur¬ 
rents  in  the  wire  objects,  and  the  eigenfunction  fields  in  the  artificial  medium.  From 
these  quantities,  the  effectives  constitutive  parameters  of  the  artificial  medium  are 
determined. 

Three  methods  are  formulated  for  determining  the  effective  constitutive  parame¬ 
ters.  The  first  method  determines  what  constitutive  parameters  a  real  medium  must 
have  to  produce  the  same  effective  wavenumber  solved  for  in  the  PMM  solution.  It 
is  shown  that  this  method  applies  only  in  certain  simple  geometries.  The  second  and 
third  methods  are  more  general,  and  are  based  on  the  assumption  that  current  and 
field  quantities  inside  the  artificial  medium  can  be  viewed  in  an  average  sense.  The 
second  method  enforces  the  constitutive  relationship  equations  using  the  current  mo¬ 
ments  and  eigenfunction  fields  averaged  over  a  lattice  cell.  The  third  method  enforces 
Maxwell’s  source-free  equations  (for  a  plane  wave)  applied  to  the  eigenfunction  fields 
averaged  over  a  lattice  cell. 

The  concept  of  artificial  media  was  first  introduced  by  Kock  [5]  in  1948  as  applied 
to  the  design  of  microwave  lenses.  The  problem  was  to  produce  a  lightweight  material 
with  suitable  index  of  refraction  at  microwave  frequencies.  Kock  presented  analysis  of 
artificial  dielectrics  consisting  of  cubic  arrays  of  spheres,  discs  and  metallic  strips.  His 
preliminary  work  did  not  include  interaction  between  the  objects,  and  thus  the  objects 
and  spacing  must  be  electrically  small  with  the  objects  small  in  terms  of  the  spacing. 
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Other  applications  include  modeling  a  plasma  as  an  artificial  dielectric  by  Rotman  [6] 
in  the  study  of  radio  wave  propagation.  Bahl  and  Gupta  [7]  designed  a  leaky-wave 
antenna  using  the  results  of  Brown  [8].  King,  Thiel  and  Park  [9]  inserted  pins  in  a 
ground  plane  to  synthesize  a  given  surface  reactance.  Sihvola  [10]  modeled  mixtures 
of  rain  and  hail  as  artificial  dielectrics  in  the  analysis  of  microwave  attenuation.  A 
recent  application  of  an  artificial  dielectric  mixture  is  in  the  microwave  welding  of 
polymers  by  Wu  and  Benatar  [11],  where  a  lossy  dielectric  of  desired  conductivity  is 
produced  by  the  proper  mixture  of  HC1  doped  polyaniline  particles  in  a  polyethylene 
host. 

Arrays  of  spherical  objects  have  been  further  treated  by  Lewin  [12]  where  he 
accounts  for  the  effects  of  plane  wave  scattering  by  nearby  metal  spheres,  and  by 
Corkum  [13]  where  he  uses  the  Clausius-Massotti  equation  [1,  Sec.  2.8.1]  suggested 
by  Kock.  Corkum ’s  analysis  accounts  the  permeability  of  such  an  array,  as  well  as 
the  permittivity  of  a  dielectric  sphere  array.  Arrays  of  thin  conducting  disks  have  also 
received  considerable  treatment.  For  the  magnetic  field  normal  to  the  disks,  Estrin 
[14]  modeled  the  current  on  an  isolated  disk  as  a  magnetic  dipole  and  solved  for  the 
permeability  in  the  case  where  the  disks  are  far  enough  apart  to  neglect  interaction. 
Estrin  [15]  also  considered  the  anisotropic  properties  of  a  3D  array  of  disks  in  his 
analysis  of  oblique  incident  waves  on  such  a  medium.  Brown  and  Jackson  [16]  include 
multipole  interaction  terms  to  provide  a  solution  accurate  for  closely  packed  disks. 
Further  analysis  of  the  metallic  strip  artificial  dielectric,  consisting  of  2D  metallic 
strips  oriented  transverse  to  both  the  direction  of  propagation  and  the  electric  field, 
is  given  by  Brown  [17]  in  which  he  formulates  a  more  accurate  theory  based  on 
transmission  line  theory.  Also,  Kolettis  and  Collin  [18]  presented  a  waveguide  modal 
analysis  for  general  directions  of  propagation  through  this  strip  media.  Experimental 
results  for  the  metallic  strip  artificial  dielectric  are  presented  by  Kolettis  and  Collin 
[18],  and  Cohn  [19]. 

Collin  [20,  Ch.  12]  presents  an  extensive  evaluation  of  artificial  dielectrics,  and 
his  work  serves  as  a  good  summary  of  much  of  the  early  work  done  in  the  area.  He 
presents  a  simple  static  Lorentz  solution  which,  only  accounts  for  dipole  interactions 
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between  objects,  and  thus,  the  objects  must  be  electrically  closely  spaced  and  small  in 
relative  to  the  spacing.  To  account  for  larger  objects,  multipole  terms  can  be  included 
in  the  expressions  for  the  fields. 

The  PMM  solution  presented  in  this  dissertation  was  first  suggested  and  utilized 
by  Blanchard  and  Newman  [21,  22].  They  analyzed  a  2D  array  of  dielectric  rods  and 
a  3D  array  of  straight  perfect  electric  conducting  dipoles.  The  current  work  included 
here  is  an  extension  of  this  preliminary  work  in  several  important  ways.  It  analyzes 
periodic  arrays  of  arbitrary  conductive  or  dielectric  wire  objects,  and  allows  for  lossy 
materials.  Also,  methods  are  formulated  to  determine  the  anisotropic  properties  of 
artificial  media,  i.e.,  the  effective  permittivity  and  permeability  tensors. 

The  theory  of  the  given  solution  to  artificial  media  is  presented  in  Chapter  2. 
This  theory  includes  the  derivation  of  the  integral  equation  for  a  plane  wave  propa¬ 
gating  in  the  artificial  medium.  The  periodic  method  of  moments  (PMM)  solution  to 
the  integral  equation  is  presented,  with  specialization  to  the  arbitrarily-shaped  thin 
wire  objects  arranged  in  the  3D  periodic  lattice  structure.  This  solution  yields  the 
complex  wavenumber,  as  well  as  the  eigenfunction  currents  and  fields,  of  the  plane 
wave  propagating  in  the  artificial  media.  A  discussion  of  how  these  quantities  are 
used  to  determine  the  effective  constitutive  parameters  of  artificial  media  is  included. 
Finally,  expressions  for  the  fields  of  a  2D  planar  array  and  a  3D  volume  array  of 
current  elements  are  derived.  These  field  expressions  provide  the  basis  for  the  PMM 
solution  and  the  eigenfunction  fields  in  the  artificial  medium. 

Chapter  3  presents  the  evaluation  of  the  MM  impedance  matrix,  and  involves  the 
computation  of  several  distinct  terms.  Chapter  4  presents  the  evaluation  of  the  aver¬ 
age  eigenfunction  electric  and  magnetic  fields  of  the  plane  wave.  The  eigenfunction 
fields  are  spoken  of  in  an  average  sense,  i.e.,  they  are  averaged  over  the  center  lattice 
cell.  Chapter  5  presents  results  obtained  from  the  solution  presented  in  this  disserta¬ 
tion.  These  results  are  intended  to  illustrate  the  methods  and  techniques,  as  well  as  to 
present  some  interesting  points  about  artificial  media.  Chapter  6  presents  the  usage 
of  the  computer  program  ADWIRS,  which  was  developed  in  FORTRAN  to  analyze 
an  artificial  medium  composed  of  a  3D  periodic  array  of  identical  arbitrarily-shaped 


thin  conductive  or  dielectric  wire  objects  arranged  in  a  homogeneous  host  medium. 
The  program  ADW1RS  yields  the  solution  of  the  plane  wave  propagating  in  the 
artificial  dielectric,  including  the  eigenfunction  currents  and  fields,  and  the  electric 
and  magnetic  dipole  moments. 

Many  of  the  mathematical  details  of  the  PMM  solution  are  given  in  appendices 
for  clarity  and  to  improve  the  readability  of  this  dissertation.  Appendix  H  presents 
a  solution  to  the  problem  of  determining  the  effective  permittivity  of  an  artifid’  di¬ 
electric  composed  of  a  3D  periodic  array  of  small  identical  dielectric  material  sp. 

The  theory  is  similar  to  that  presented  in  Chapter  2,  but  it  is  not  as  general  and 
complete.  However,  the  implementation  of  the  method  is  quite  different,  due  to  the 
scattering  objects  being  spheres  instead  of  wire  objects.  Appendix  H  is  intended 
to  be  a  self-contained  document,  with  few  references  to  material  elsewhere  in  this 
dissertation.  -/ 
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Chapter  2 
Theory 


This  chapter  presents  the  integral  equation  and  PMM  solution  for  a  plane  wave  prop¬ 
agating  through  an  artificial  medium  composed  of  thin  conductive  or  dielectric  wire 
objects  arranged  in  a  periodic  lattice.  This  solution  yields  the  complex  wavenum¬ 
ber  of  the  plane  wave  that  can  propagate  without  excitation  through  the  artificial 
medium,  i.e.,  the  eigenfunction  solution  for  the  artificial  medium.  The  solution  also 
yields  the  shape  of  the  eigenfunction  current*  in  the  scattering  objects,  from  which 
the  eigenfunction  fields  can  be  determined.  Throughout  this  paper,  all  fields  and 
currents  are  assumed  to  be  time  harmonic  with  the  e**  time  dependence  suppressed. 

Two  methods  are  presented  for  determining  the  effective  permittivity  and  perme¬ 
ability  tensors  of  the  artificial  medium.  The  polarisation  method  is  based  on  enforcing 
the  constitutive  relationship  equations  in  an  average  sense  in  the  artificial  medium, 
and  uses  the  average  eigenfunction  currents  and  fields  per  unit  volume.  The  Maxwell’s 
equations  method  is  based  on  the  assumption  that  Maxwell’s  equations  for  the  plane 
wave  are  satisfied  in  an  average  sense  in  the  artificial  medium.  This  method  uses  the 
average  eigenfunction  fields  per  unit  volume  and  the  complex  vector  wavenumber  of 
the  plane  wave. 

As  shown  in  Figure  2.1,  the  geometry  of  the  artificial  medium  consists  of  a  3D  triple 
periodic  array  of  wire  objects  located  in  a  homogeneous  and  isotropic  host  medium. 
The  homogeneous  host  medium  has  constitutive  parameters  (p o,«o)»  wavelength  Ao, 
wavenumber  ho,  and  is  not  necessarily  free  space  and  may  be  lossy.  The  wire  objects 
(shown  as  V-dipoles  in  the  Figure  2.1)  may  be  composed  of  an  arbitrary  conductive 
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or  dielectric  material.  The  wires  have  radius  a,  constitutive  parameters  (/io,ei),  and 
wavenumber  hi.  The  complex  permittivity  of  the  wire  objects,  Cj,  is  given  by 

«i  =  cir  eo  -  j  —  =  €lr  co  (1  -  j  tan  Si)  (2.1) 

w 

where  ox  is  the  conductivity  of  the  wire  object  (in  ft-1/  meter)  and  tanfi  is  the  loss 
tangent  of  the  wire  objects.  The  wire  objects  are  arranged  in  a  lattice  cell  structure 
with  spadngs  dy,d„  and  d*,  in  the  fi,$  and  w  directions,  respectively.  This  lattice 
cell  structure  need  not  be  rectangular,  and  its  unit  vectors  are 

6  =  *, 

v  =  v,x  +  t/*y  and 

w  =  to*  ±  +  wv  f  +  wa  i. 

The  elements  are  referenced  by' the  index  Q  =  (•*,  jv,K>)  where  — oo  <  (ttt,jv,h»)  < 
oo.  The  reference  or  center  wire  object  is  centered  at  the  origin  and  is  indexed  by 
Q  =  0  =  (0, 0, 0).  Also,  let  AQ  =  »Mdyfi + j0d„V  -f  Kid^r  be  the  position  vector  from 
the  origin  to  the  center  of  lattice  cell  Q. 

Following  the  general  methods  presented  by  Blanchard  and  Newman  [21,  22],  the 
effective  parameters  of  the  artificial  medium  are  determined  by  first  assuming  that 
a  plane  wave  of  unknown  wavenumber  is  propagating  through  the  medium.  If  the 
plane  wave  is  propagating  in  the  U*  direction,  assumed  to  be  known,  then  it  will  be 
of  the  form 

e-*«R,  (2.2) 

where  k «  =  fee  fi*  =  k^x  +  k^y  +  ke,  %  is  the  unknown  complex  vector  wavenumber 
(or  wave-vector) j  and  R  =  xx  -(-  yf  +  xi  is  the  position  vector.  For  a  given  direction 
fit,  it  is  desired  to  find  K  such  that  this  plane  wave  satisfies  Maxwell’s  source  free 
equations  and  all  of  the  boundary  conditions  in  the  artificial  medium,  corresponding 
to  the  normal  mode  of  propagation  for  the  artificial  medium,  i.e.,  the  eigenfunction 
solution  for  the  artificial  medium.  Once  the  wavenumber  &e  is  known,  then  the  shape 
of  the  eigenfunction  currents  in  the  material  wire  objects  can  be  determined.  Rom 
these  eigenfunction  currents,  the  eigenfunction  fields  can  then  be  found. 
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Wire  Object 

Q  =  (iutjvikw) 

(Vo^i) 


v=v*x+vyy 


Background 

Medium 

( Mo»eo) 


w=wxx+wyy+wzz 


Figure  2.1:  The  3D  volume  lattice  geometry. 


For  plane  wave  propagation  in  a  given  direction  through  a  homogeneous  aniso¬ 
tropic  medium,  it  is  known  that  there  are  two  distinct  eigenfunction  modes  of  prop¬ 
agation  [23,  Sec.  4.25],  [24,  Sec.  14.2.2].  Each  mode  of  propagation  corresponds  to 
a  distinct  complex  wavenumber  and  polarisation  of  the  plane  wave.  Also,  for  homo¬ 
geneous  anisotropic  media,  the  nine  elements  of  the  permittivity  and  permeability 
tensors  are  independent  of  the  direction  of  propagation  through  the  medium  [23,  Ch. 
4],  [24,  Sec.  14.1]. 

It  will  be  shown  that  the  theory  of  two  distinct  modes  of  plane  wave  propaga¬ 
tion  extends  to  plane  wave  propagation  through  artificial  media.  Thus,  for  a  given 
direction  of  propagation  u*  through  an  artificial  medium,  there  are  two  values  of 
A,  corresponding  to  the  two  eigenfunction  solutions  for  the  artificial  medium.  Each 
value  of  he  corresponds  to  a  distinct  polarisation  of  the  plane  wave,  and  to  distinct 
eigenfunction  currents  on  the  wire  objects,  and  hence,  to  distinct  eigenfunction  fields 
in  the  artificial  medium.  Also,  it  will  be  shown  that  the  permittivity  and  permeability 
tensors  typically  will  have  only  a  slight  variation  on  the  direction  of  propagation  in  an 
artificial  medium.  One  possible  exception  is  when  the  eigenfunction  current  shapes 
are  a  strong  function  of  the  direction  of  propagation. 

2.1  The  Integral  Equation 

In  formulating  the  integral  equation  for  the  artificial  medium,  as  shown  in  Figure  2.2, 
the  volume  equivalence  theorem  is  used  to  replace  the  wire  objects  by  the  host  medium 
and  the  equivalent  electric  volume  polarisation  currents  [1,  Sec.  7.7] 

J  =  yu>(ci  —  eoJE*,  (2*3) 

where  E*  is  the  total  electric  field  inside  the  wire  objects.  Since  the  permeabilities  of 
the  host  medium  and  the  wire  objects  are  identical,  there  are  no  magnetic  currents  in 
the  wire  objects.  In  the  limiting  case  where  the  wires  become  perfectly  conducting, 
the  volume  current  J  approaches  a  surface  current  on  the  wire  surfaces.  The  current 
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J  exists  in  (or  on)  each  and  every  wire  object,  and  is  written  as 

J  =  (2.4) 

Q 

where  jQ  is  the  current  in  wire  object  Q,  and  the  summation  is  over  all  values  of 
Q,  i.e.,  —oo  <  (iu,jv,K>)  <  oo.  Since  we  seek  a  solution  to  Maxwell's  source  free 
equations,  there  are  no  impressed  currents,  and  thus  E*  is  the  electric  field  of  J 
radiating  in  the  homogeneous  host  medium.  Equation  (2.3)  can  be  rearranged  as  the 
homogeneous  equation 

.  J 

—  E*  H-  *7-7 - r  =  0  in  each  wire  object  Q,  (2.5) 

-  eo) 

and  is  to  be  solved  for  the  complex  wavenumber  fee,  and  the  current  in  the  center  or 
Q  =  0  wire  object,  by  the  PMM. 

Due  to  the  periodic  nature  of  the  array  of  wire  objects,  and  of  the  plane  wave  of 
Equation  (2.2),  the  current  is  identical  in  each  wire  object  except  for  an  amplitude 
and  phase  change  corresponding  to  the  value  of  the  plane  wave  at  the  center  of  the 
lattice  cell.  In  other  words,  the  current  in  wire  object  Q  differs  from  the  current  in 
the  center  wire  object  by  the  complex  multiplier 

CQ  =  e~*«-AQ.  (2.6) 

As  a  result,  the  only  unknowns  are  K  and  the  current  in  the  center  wire  object. 

2.2  PMM  Solution  of  the  Integral  Equation 

Equation  (2.5)  will  be  solved  by  the  periodic  moment  method  (PMM)  [2,  3,  4].  The 
first  step  in  the  PMM  solution  is  to  expand  the  unknown  current  J  as 

j -£**•»  (2-7) 

q  q  n=i 

where  the  J|?  are  N  linearly  independent  expansion  functions  for  the  current  in  wire 
object  Q,  and  the  I„  are  N  unknown  expansion  coefficients,  with  n  =  1,2, ... ,  JV. 
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Equivalent 


Figure  2.2:  The  equivalent  3D  volume  lattice  geometry  after  application  of  the  volume 
equivalence  theorem. 


Note  that  for  a  given  n,  all  the  Jjp  are  identical  in  shape,  with  the  only  difference 
bong  that  a  given  J?  is  defined  only  on  wire  object  Q,  i.e., 

•>!(*)  =  J?(R  -  AQ).  (2.8) 

Due  to  the  periodic  nature  of  the  problem,  it  is  only  necessary  to  enforce  Equa¬ 
tion  (2.5)  over  V° ,  the  volume  of  the  center  wire  object. 

Next,  define  N  linearly  independent  weighting  functions  in  the  center  wire  object, 
denoted  as  Wm  with  m  =  1,2, ...,  N.  Substituting  J  of  Equation  (2.7)  into  Equa¬ 
tion  (2.5),  and  taking  the  inner  product  of  the  result  with  the  N  weighting  functions, 
reduces  Equation  (2.5)  to  an  order  N  matrix  equation  which  can  be  written  as 

[Z  +  AZ]J  =  [Z]J  =  0.  (2.9) 

Here,  [Z]  —  [Z  - f  A Z]  is  the  6rder  N  impedance  matrix  and  I  is  the  length  N  solution 
vector  containing  the  /„  expansion  coefficients  of  Equation  (2.7).  The  impedance 
matrix  elements  are  given  by  (m,n  =  1,2,...,  N) 

=  £  C<*z2»  =  -£<?“/  E?  ■  W„  iv  (2.10) 

Q  Q  m 

AZ„„  =  — i - rf  (2.11) 

je>(e  i-€o)Jvm 

where  is  the  electric  field  of  J^,  the  nth  expansion  function  in  wire  object  Q, 
radiating  in  the  homogeneous  host  medium.  The  integration  is  over  Vm,  the  volume 
of  weighting  function  m.  The  AZmn  terms  do  not  depend  on  K,  but  the  Zmn  terms 
do  through  their  dependence  upon  C^.  Note  that  the  A Zmn  terms  are  non-sero  only 
when  the  expansion  functions  J?  are  in  the  center  Q  =  0  wire  object.  Also,  the 
AZmn  terms  vanish  when  the  wires  are  perfectly  conducting,  since  |ci|  — ►  0.  The 
impedance  matrix  terms,  Zmn(ke)  and  A Zmn,  are  evaluated  in  Chapter  3,  for  the 
piecewise  sinusoidal  (PWS)  material  wire  expansion  and  weighting  functions  defined 
in  Section  2.2.1. 

The  homogeneous  matrix  Equation  (2.9)  will  have  a  non-trivial  solution  only  if 
the  determinant  of  the  impedance  matrix  is  zero.  Thus,  K  is  found  by  solving  the 
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fundamental  equation 

I  Z(k.)  +  AZ|  =  \Z\  =  0,  (2.12) 

usually  on  an  iterative  basis. 

2.2.1  The  MM  Expansion  and  Weighting  Functions  for  Ma¬ 
terial  Wires 

In  this  section,  MM  expansion  and  weighting  functions  suitable  for  thin  conductive 
or  dielectric  wires  are  defined.  The  expansion  functions  used  are  those  employed  by 
Newman  [25],  which  incorporate  known  behavioral  variations  of  the  thin  material 
wires.  As  shown  in  Figure  2.3,  a  dipole  expansion  function  n  consists  of  two  monopole 
current  segments.  Monopole  nj  is  oriented  along  the  znj-axis,  with  zero  current 
at  znj  =  0  rising  to  unity  current  at  Znj  =  dnj.  The  polarity  of  the  current,  on  a 
monopole  basis,  always  flows  from  the  zero  current  end  toward  the  unity  current  end. 
However,  dipole  expansion  function  n  has  polarity  such  that  current  flows  from  the 
j  =  1  monopole  to  the  j  =  2  monopole.  Current  is  continuous,  and  of  value  unity, 
across  the  terminals  of  the  dipole.  The  terminals  of  the  dipole  is  where  monopole 
1  intersects  with  monopole  2,  i.e.,  the  point  znj  —  d„j.  In  this  manner,  expansion 
function  n  is  written  as  (n  =  1,2,...,  AT) 

Jn  =  > znl )  *n2)  (2.13) 

where,  for  j  =  1,2 

Jnj(Pnj,  Znj)  =  ^njiPnjt  znj)  +  pnj  Jnj(pnj>  ^j)  (2*14) 

with  p„j  and  z„j  bang  the  local  radial  and  axial  coordinates  for  monopole  segment 
nj. 

The  expansion  function  contains  an  axial  and  a  radial  component  corresponding 
to  the  transverse  magnetic  (TM)  to  z  fields  inside  the  thin  material  wire  [25].  These 
components  are  written  as 


Jt. IjiPnjtZnj)  —  C  Jo(kppnj)  Fnj[hz  Znj)  and 


(2.15) 


Jnj(pnjyznj)  ~  r  *^1  ( kp Pnj )  Fnj ( kM  Z„j ) 

Kp 

where  C  is  a  normalization  constant,  J,  is  a  Bessel  function  of  order  t  (t  =  0, 1),  the 
prime  '  denotes  differentiation  with  respect  to  z„j,  and  from  the  wave  equation, 

k]  +  k]  =  k\. 


The  normalization  constant  C  is  chosen  such  that  znj)  has  unit  terminal 

current,  i.e., 

C=il  L  H  =  2Ta  Ji{kp a) ’  (2’l7) 

Note  that  the  radial  component  J^j  is  dependent  upon  the  axial  component  J*j,  and 
thus  there  is  only  one  unknown  current.  The  function  Fnj(kx  znj)  defines  the  axial 
variation  of  the  expansion  functions,  as  explained  next. 

As  shown  in  Figure  2.3,  dipole  expansion  functions  are  always  zero  at  the  endpoints 
but  rise  to  unity  at  the  terminals.  Since  the  current  vanishes  at  the  endpoints  of 
a  perfectly  conducting  wire,  these  expansion  functions  are  used  for  modeling  the 
current  on  perfectly  conducting  wire  objects.  They  are  also  used  for  modeling  the 
current  away  from  wire  endpoints  on  imperfect  conductive  or  dielectric  wires  objects. 
The  axial  variation  of  a  dipole  expansion  function  consists  of  two  monopole  axial 
variations.  Choosing  kx  =  ko,  the  monopole  axial  variation  is 


if  0  <  znj  <  dnj 
0  otherwise. 


«in(*o  *ni) 

«in(*o  dn/) 


(2.18) 


For  imperfect  conductive  or  dielectric  wire  objects,  or  when  wire  objects  physically 
touch  across  adjacent  cells,  the  current  will  not  necessarily  vanish  at  the  wire  object 
endpoints.  To  model  such  a  current,  monopole  expansion  functions  are  used.  A 
monopole  expansion  function  consists  of  only  the  j  =  1  current  segment.  In  analogy 
to  Equation  (2.13),  monopole  expansion  function  n  is  written  as 


dn  —  Jnl(pnl  j  2fil  )j 


(2.19) 
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Zn2=dn2 


Figure  2.3:  A  typical  dipole  MM  expansion  function. 


and  the  axial  and  radial  variation!  remain  the  same  as  for  the  dipole  expansion 
functions. 

The  weighting  functions  contain  only  an  axial  component  and  have  axial  varia¬ 
tion  the  same  as  the  expansion  functions.  However,  the  radial  variation  is  constant. 
Therefore,  a  dipole  weighting  function  is  defined  as 

=  I  ^ml)  ®m2  ^m2 (h«  Zm 2)  ]  (2.20) 

and  a  monopole  weighting  function  is  defined  as 

Wm  =  -^amlFml(*,zml)  (2.21) 

for  m  =  1,2, ...,  JV. 


2.2.2  Evaluation  of  the  Eigenfunction  Currents 


Assume  that  a  value  of  kg  has  been  found  such  that  the  fundamental  Equation  (2.12) 
is  satisfied.  Now  it  is  desired  to  determine  the  eigenfunction  currents  in  the  wire 
objects.  The  MM  matrix  equation  for  the  eigenfunction  current  in  the  center  wire 
object  can  be  written  as 


Zu  Z13  •  •  •  Z\n 

> 

0 

Z22  •  •  •  Z2N 

h 

0 

*  •  •  • 

•  •  •  • 

•  •  •  • 

Zm  Zn2  •  •  •  Znn 

•  •  •  .  ^ 

• 

• 

• 

0 

(2.22) 


Recall  that  the  determinant  of  the  impedance  matrix  is  zero,  and  thus  this  system  of 
equations  cannot  be  solved  for  the  current  coefficients  I  in  the  usual  manner.  However, 
the  eigenfunction  currents  can  be  determined,  to  within  a  constant,  by  setting  an 
arbitrary  non-zero  element  of  I  to  unity,  and  then  reducing  Equation  (2.22)  to  an 
order  N  —  1  system  of  equations  for  the  remaining  N  —  1  current  coefficients.  For 
example,  setting  the  lafet  coefficient  In  =  1,  Equation  (2.22)  reduces  to 


Z\\  Z\2  •  •  •  Z\t  N-\ 

/1 

—Z\N 

Z21  Z22  •  •  •  Z2,N-\ 

h 

~Z\N 

•  •  •  • 

•  •  •  • 

•  •  •  • 

_  Zn-  1,1  Zn-  1,2  •  •  •  Zn-\,n-\ 

• 

• 

• 

In-  1  _ 

* 

H 

1 

...  > 

1 

(2.23) 
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which  can  be  solved  using  standard  techniques.  This  method  yields  the  eigenfunction 
currents  on  the  wire  objects  to  within  a  constant. 

2.3  Determination  of  the  Effective  Constitutive 
Parameters 

This  section  discusses  the  determination  of  the  effective  permittivity  and  permeability 
tensors,  denoted  by  (ee,/ie),  for  an  anisotropic  artificial  medium.  A  discussion  of  the 
roots  he  to  Equation  (2.12)  is  included,  and  different  approaches  to  evaluating  (ce,jie) 
are  presented. 

Much  of  the  theory  contained  in  this  section  follows  from  anisotropic  media  theory 
as  applied  to  crystal  optics  [23,  Ch.  4],  [24,  Ch.  14].  This  theory  seems  to  apply 
to  the  macroscopic  model  of  artificial  media,  provided  the  arbitrary  material  wire 
objects  are  not  too  electrically  large  or  spaced  too  far  apart  and  the  eigenfunction 
current  shape  is  not  a  strong  function  of  the  direction  of  propagation. 

2.3.1  Discussion  of  the  Roots  ke  and  Polarization 

The  first  step  in  determining  the  equivalent  permittivity  and  permeability  of  the  arti¬ 
ficial  medium  is  to  determine  the  roots,  he,  to  Equation  (2.12).  For  a  given  direction 
of  propagation  through  an  anisotropic  artificial  medium  composed  of  arbitrary  scat¬ 
tering  objects,  there  are  two  fundamental  roots,  he,  to  Equation  (2.12)  [23,  Sec.  4.25], 
[24,  Sec.  14.2.2].  Each  root  corresponds  to  a  distinct  polarization  of  the  plane  wave 
in  the  artificial  medium,  and  its  corresponding  eigenfunction  currents  and  fields. 

For  special  geometries,  the  roots  may  be  repeated  roots  or  degenerate  to  the  host 
medium  wavenumber  he  =  ho.  For  example,  repeated  roots  (two  roots  with  the 
same  numerical  value  for  he)  will  occur  for  propagation  normal  to  a  symmetric  wire 
cross  with  equal  length  vertical  and  horizontal  members.  One  root  corresponds  to 
vertical  polarization,  and  the  other  root  corresponds  to  horizontal  polarization.  If  the 
scattering  object  is  of  2D  or  of  3D  extent,  the  two  roots  will  be  different  from  the  host 
medium  wavenumber.  If  the  scattering  object  is  of  ID  extent  (i.e.  a  linear  dipole), 
then  there  will  be  one  root  he  =  h o,  corresponding  to  a  plane  wave  with  polarization 
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perpendicular  to  the  dipole,  and  one  root  he  ^  Jfeo,  corresponding  to  polarisation 
parallel  to  the  dipole. 

In  a  real  homogeneous  anisotropic  medium  the  elements  of  the  tensor  constitutive 
parameters  are  independent  of  the  direction  of  propagation  [23,  Ch.  4],  [24,  Sec.  14.1]. 
However,  in  an  artificial  medium,  these  elements  may  depend  on  the  direction  of  prop¬ 
agation.  If  the  material  wire  objects  are  not  too  electrically  large  or  spaced  too  far 
apart,  and  the  eigenfunction  current  shape  is  not  strongly  dependent  on  the  direction 
of  propagation,  then  any  direction-dependent  variation  of  the  tensor  constitutive  pa¬ 
rameter  elements  should  be  small  or  negligible.  However,  if  the  eigenfunction  current 
shape  is  a  strong  function  of  the  direction  of  propagation,  then  the  tensor  constitutive 
parameter  elements  may  also  vary  with  the  direction  of  propagation. 

Note  that  in  developing  Equation  (2.12),  the  polarization  of  the  plane  wave  was 
not  specified.  Thus,  when  a  foot  of  Equation  (2.12)  is  found,  the  polarization  is  un¬ 
known.  To  find  the  polarization  one  can  first  compute  the  corresponding  eigenfunc¬ 
tion  currents  using  Equation  (2.23),  and  then  the  eigenfunction  electric  and  magnetic 
fields,  as  presented  in  Chapter  4.  The  polarization  of  the  electric  field  is  taken  as  the 
polarization  corresponding  to  the  chosen  value  of  k-  Once  the  eigenfunction  fields 
are  known,  the  characteristic  impedance  of  the  artificial  medium,  denoted  T)e,  can  be 
evaluated  as  the  ratio  of  the  electric  to  magneto  eigenfunction  fields  tangential  to 
the  assumed  direction  of  propagation.  If  the  host  medium  and  the  scattering  objects 
are  lossless,  then  k  and  i\e  will  be  positive  real  numbers.  However,  if  either  the  host 
medium  or  the  scattering  objects  have  loss,  then  k  will  be  a  complex  number  in  the 
fourth  quadrant,  and  rfe  will  be  a  complex  number  in  the  sector  ±45°  of  the  positive 
real  axis  [26,  Sec.  2-3]. 

Since  the  scattering  objects  in  an  artificial  medium  are  typically  electrically  small, 
one  is  usually  interested  in  the  eigenfunction  modes  with  the  smallest  k>  However,  it 
is  important  to  note  that  larger  roots,  with  larger  he,  may  exist.  One  manner  in  which 
higher  order  roots  do  exist  is  through  the  periodic  nature  of  the  complex  multiplier 
CQ  of  Equation  (2.6).  For  example,  in  a  problem  where  propagation  is  in  the  u*  =  i 
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direction,  simplifies  to 


CQ  =  e~jk,(kwdm  w,)' 


(2.J4) 


Recall  that  is  the  only  term  in  the  fundamental  root  Equation  (2.12)  that  depends 
on  kf.  Let  K  =  K o  be  the  smallest  root  to  Equation  (2.12),  where  Re(fceo)  >  0.  It  can 
be  seen  by  inspection  of  Equation  (2.24)  that  if  fceo  is  a  root,  then  so  will  be  values 
of  fee  where 

K  =  Ko±^~  for  p=  1,2,3,...  (2.25) 

d^w. 

Thus,  for  some  geometries,  there  exists  a  periodic  occurrence  of  roots  to  Equa¬ 
tion  (2.12).  However,  for  the  material  presented  here,  the  lowest  order  root  is  the 
only  root  of  interest,  i.e.,  K  =  fceo  is  assumed.  It  should  be  noted  that  higher  order 
roots  may  exist  that  correspond  to  higher  order  modal  eigenfunction  solutions  in  the 
artificial  medium,  however,  no  such  roots  have  been  found  to  this  date. 


2.3.2  Polarization  Method 


Consider  the  determination  of  the  dyadic  effective  permittivity  and  permeability 
(ce,fte)  for  an  anisotropic  artificial  medium.  It  is  assumed  that  for  a  given  geom¬ 
etry,  the  two  roots  Afe,  their  corresponding  eigenfunction  currents  J°  on  the  center 
element,  and  their  eigenfunction  fields  averaged  over  the  volume  of  the  center  cell 
(E°,H°),  have  all  been  determined.  In  the  limit  as  ke  -*  ko,  the  eigenfunction  cur¬ 
rents  vanish  and  the  eigenfunction  fields  become  identical  to  a  plane  wave  in  the  host 
medium. 

The  electric  and  magnetic  dipole  moment  per  unit  volume  in  the  center  cell  can 
now  be  computed  as  [20,  Ch.  12.5] 

P'  =  ]^,h  (2-26) 

M° =  sL /vR  *  3' dv  =  *”H*  <2-2r> 

where  F°  represents  the  region  of  the  center  cell  of  volume  Av,  and  x*  and  xm 
are  the  dimensionless  symmetric  [23,  Sec.  4.24],  [24,  Sec.  14.1]  effective  electric 
and  magnetic  susceptibility  tensors,  respectively,  for  the  artificial  medium.  Equation 
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(2.27)  shows  that  even  perfectly  conducting  or  pure  dielectric  scattering  objects  can 
have  a  magnetic  moment,  and  thus  an  effective  permeability  different  from  the  host 
medium.  Assuming  that  the  usual  constitutive  relationships,  which  are  valid  point  by 
point  in  a  real  medium,  hold  in  an  average  sense  in  an  artificial  medium,  the  average 
electric  and  magnetic  flux  densities  in  the  center  cell  are  given  by 

D®  =  £oE°  +  P®  =  eoE®  +  «o*e  •  E®  =  ce  •  E®  (2.28) 

B®  =  /io(H®  +  M°)  =  /*o(H®  +  xm  •  H®)  =  Ae  *  H®.  (2.29) 

From  Equations  (2.28)  and  (2.29),  the  dyadic  effective  permittivity  and  permeability 
are  given  by 

=  (/  +  **)« 0  (2.30) 

Ae  =  (*  +  Xm)*i  (2.31) 

J 

where  I  is  the  unit  dyad. 

Explicitly  showing  Equation  (2.26)  relating  the  effective  electric  susceptibility  to 
the  average  electric  field  and  the  electric  dipole  moment  per  unit  volume  in  the  center 
cell, 

'*.*,*.][<] 

“•444  sj  -  p!  •  (2-32) 

xl*  xL  -®?  p‘ 

Equation  (2.32)  is  equivalent  to  three  equations  in  the  nine  components  of  %*,  and  is 
the  result  of  one  of  the  two  roots  of  Equation  (2.12).  The  other  root  will  produce  a 
dyadic  equation,  similar  to  Equation  (2.32)  with  the  same  x®,  but  with  different  E® 
and  P®.  The  two  dyadic  equations,  along  with  the  condition  that  xe  symmetric, 
can  now  be  solved  for  the  nine  components  of  \e .  Once  \e  »  known,  then  ie  is 
determined  simply  from  Equation  (2.30).  The  determination  of  Ae  i®  parallel  to  that 
presented  for  ee,  but  uses  the  dyadic  Equation  (2.27). 

2.3.3  Maxwell’s  Equations  Method 

Another  method  of  determining  the  tensor  constitutive  parameters  of  the  artificial 
medium  follows  directly  from  Maxwell’s  equations.  For  this  method,  it  is  assumed 
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that  Maxwell’s  equations  apply  to  the  average  fields  in  the  center  cell  of  the  artificial 
medium.  Applying  the  well-known  source-free  Maxwell’s  equations  to  the  plane  wave 
fields  of  the  form  assumed  in  Equation  (2.2),  it  is  obtained  that  [27,  Sec.  2.3] 


VxH  =  -ju>D  =>  It,  x  H®  =  -wD1  =  -we,  •  E#  (2.33) 
V  x  E  =  jw B  ^  k,  x  E®  =  wB®  =  u>fte  •  H®  (2.34) 

V  •  D  =  0  =►  ke  •  D®  =  0,  (2.35) 

VB  =  0  =*>  lceB®  =  0  (2.36) 


where  the  tensor  constitutive  relationships  of  Equations  (2.28)  and  (2.29)  have  been 
used. 

Explicitly  showing  Equation  (2.33)  relating  the  effective  wave-vector  and  permit¬ 
tivity  tensor  to  the  average  electric  and  magnetic  fields  per  unit  volume  in  the  center 
cell,  ' 

[4,  4,  4,1  [i? 

=-"  4,  4,  4.  E;  ■  (2-37) 

K,Hi-K,a:\  [4,  4,  <.][lf_ 

Equation  (2.37)  is  equivalent  to  three  equations  in  the  nine  components  of  ee,  and  is 
the  result  of  one  of  the  two  roots  of  Equation  (2.12).  The  other  root  will  produce  a 
dyadic  equation,  similar  to  Equation  (2.37)  with  the  same  ee,  but  with  different  k«  and 
(E°,H°).  The  two  dyadic  equations,  along  with  the  condition  that  £e  «  symmetric, 
can  now  be  solved  for  the  nine  components  of  ee.  The  determination  of  fie  is  parallel 
to  that  presented  for  !e,  but  uses  the  dyadic  Equation  (2.34). 

2.3.4  Uniaxial  Artificial  Media 

For  certain  simple  anisotropic  media,  termed  uniaxial  media,  the  off  diagonal  compo¬ 
nents  of  x*  (and  xm)  ***  negligible,  and  the  diagonal  components  are  related  to  E* 
and  P®  by 

X£  =  ^|p  »  =  *.?,*•  (2.38) 

Equation  (2.38)  can  be  used  to  find  the  ii  component  of  jjf  provided  that  E*  is  non 
aero.  If  for  a  given  direction  of  propagation  and  root,  E?  ^  0  for  t  =  z,y,z,  then 
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Equation  (2.38)  can  be  used  to  find  all  three  diagonal  components  of  x*.  If  for  a 
given  direction  of  propagation  and  root  £?  =  0,  then  P*  =  0,  and  Equation  (2.38) 
is  indeterminate.  In  this  case  x?  can  be  determined  from  the  other  root,  or  by  a 
different  direction  of  propagation. 

Consider  propagation  along  of  one  of  the  three  principle  axes  in  a  uniaxial  media. 
The  two  roots  will  correspond  to  polarisations  in  the  directions  of  the  two  principle 
axes  transverse  to  the  direction  of  propagation.  For  this  special  case,  the  effective 
permittivity  and  permeability  can  be  determined  from  the  root  he  and  characteristic 
impedance  qe.  It  is  assumed  that  ke  and  rje  have  been  determined  for  a  given  frequency, 
polarisation  and  direction  of  propagation  along  one  of  the  principle  axes.  Since  K 
and  Tje  are  related  to  /te  and  e«  by 


K  t  “d  Vt  = 


(2.39) 


then 


fce  j  tyefce 

e*  = -  and  /x*  = - . 

Wife  & 


(2.40) 


If  the  magnetic  moment  of  the  scattering  objects  is  negligible,  then  fie  =  Ho,  and  e* 
is  given  by 

(2.41) 


Jfc» 

Ce  =  . 

W3Ho 


The  advantage  of  using  Equation  (2.41)  for  pure  artificial  dielectrics  (i.e.,  fie  =  fio)  is 
that  one  need  not  compute  ife. 


2  A  Fields  of  a  2D  Planar  Array  of  Current  Ele¬ 
ments 

This  section  evaluates  the  exact  electric  and  magnetic  fields  of  a  2D  planar  array  of 
linear  current  elements.  The  theory  included  here  has  been  presented  before  [28,  29, 
30],  however,  it  is  repeated  here  since  it  is  of  great  importance  in  implementing  the 
PMM  solution,  and  to  employ  notation  specific  to  the  artificial  dielectric  geometry. 

Figure  2.4  shows  a  2D  planar  array  of  infinitesimal  current  elements  arranged  on  a 
parallelogram  grid.  The  host  medium  is  homogeneous  and  isotropic  with  constitutive 
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parameters  (fio,  eo).  Note  that  the  host  medium  is  not  necessarily  free  space,  and  in 
general  is  a  lossy  medium.  The  unit  vectors  of  the  grid  are  d  =  x  and  v  =  vxx  -f  vwy, 
and  the  spadngs  of  the  grid  are  d*  in  the  u  direction  and  d*  in  the  v  direction.  The 
elements  are  referenced  by  the  indices  (iu,jv)  for  — oo  <  (s«,y«)  <  oo  and  element 
(iu,jv)  is  located  at  R,j  =  iu  d*  u + jv  dy  v.  All  elements  are  polarised  in  the  arbitrary 
direction  a,  which  may  have  a  component  out  of  the  xy  plane  (or  the  plane.) 
Corresponding  to  the  value  of  the  plane  wave  with  wave-vector  k«  propagating  across 
the  array,  the  incremental  current  moment  on  element  (*„,  jv)  is 


kldl'  =  kldl'  e-X*  *+*  where 

k«  =  kesx  +  keyy  +  keti, 
di  =  dy  keX,  and 

dj  =  MketVs  +  kvVy). 

The  entire  array  produces  some  incremental  electric  vector  potential  dA,  which  is 
the  the  summation  of  each  individual  element's  contribution.  If  the  field  point  is  at 
R  =  xx  +  yy  +  si,  then  the  contribution  of  element  ( iu,jv )  to  dA  is 

/* o 


dAtj  —k^-Idl  e~X,m *+» *>  S.  . 

u  4*  |R  -  Rij\ 

Summing  the  contributions  of  every  element  in  the  array,  it  is  obtained  that 


u  •  « 


-i(i»  di)  V" 


(2.42) 


(2.43) 


..  |It  -  • 

In  Equation  (2.43)  it  is  implicit  that  the  summations  are  over  —  oo  <  (*u,j„)  <  oo. 
The  inner  summation  of  Equation  (2.43)  can  be  rewritten  as 


_  ...  <*«)*+«* 
Y'  <*.)  _ _ _  v  *) 

^  |R  -  R*l  t"  v((*' +  ®1 


(2.44) 

//  — ./  -  J  12  I  —JL 

•  «  I  •« 

where  a'  —  x  —  jv  vx  dv  and  a3  =  (v  —  i.vvd.)3  +  23.  Next,  the  Polison  nun  formula 
[31] 

£e*^‘i-(mm,)  =  -£/(<  +  »-)  (2.45) 

m  w0  n  '  W0  / 

along  with  the  Fourier  transform  pair 

_-ifcoVa,+("-"i)*  >i(  ...  /  / - \ 

f(“} ■  ^ ~ '<«> ~  !Tg°  (a^) ’  (2'w) 
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where  JJ^(*)  is  the  second  type  Henkel  function  of  order  sero,  is  applied  to  Equa- 
tion  (2.44).  With  the  variables  of  transformation  defined  as 

m  — *  »„ 

»  -*  n, 

«o  -»  d* 

v i 

i  — ►  —he* 

it  is  obtained  that 

V  e  _  *  y  ->*(*«.+*»  £)  J(j.«.A.)(fc*<,+n«  35)  V 

t  iR-Riii  "id-tr 

x  H<‘>  afa-^  +  ^y  .  (2.47) 

i 

All  the  summation  indices  of  Equations  (2.45)  and  (2.47)  are  of  the  range  — oo  < 
(m,n,nu)  <  oo 

The  electric  vector  potential  can  be  now  be  written  as 

iK  =  *e2^£e-*>*<fr+"-*>x 

ft* 

x  £  IkorJiy-j^f  +  zA  ,  (2.48) 

i»  L  J 

where  rp  =  ^1  —  (■%£■  -f  The  Poisson  sum  formula  of  Equation  (2.45)  is 

applied  a  second  time  using  the  Fourier  transform  pair 

F(u)  =  (hoT^Z*  +  (u,  -  a*)*)  <=►  f(t)  = 

and  the  transformation  variables  defined  as 

m  -♦  jv 
n  — ►  n» 

Mo  -»  Vyd* 

«2  -*  V 

t  -*0  (fe1 -«.“&)  • 
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Carrying  out  this  process  on  the  inner  summation  of  Equation  (2.48),  it  is  obtained 
that 


dA  =  & . £qI*L.  y +"-»&) y 

2jk0duvvdv  X  T? 


^  -  (fr  + n“ft)2  ~(i£-  + 

It  is  worthwhile  to  write  Equation  (2.50)  in  the  compact  form 


(2.50) 


J4  .  Holdl'  *^  •-*■**  , 

dA  =  a;m — t  2-  X, -  for  2 

2jk0duvydv^^t  r* 


where 


(2.51) 


R  =  *£  +  If?  +  *i, 
r±  =  rxx  +  r¥y  ±  r,s, 
r,  = 

r»  =  if  ~  n“i£fcfr + n°  *o^i. » 

r,  =  ^/l  —  r*  —  r*  such  that  Im(r, )  <  0. 

In  this  form,  the  electric  field  of  the  2D  array  of  infinitesimal  current  elements  is 

expressed  as  a  double  spectral  summation  of  plane  waves  propagating  in  the  spectral 

direction  r+  for  positive  z  and  f_  for  negative  z.  Note  that  in  the  square  root  defining 

r„  the  root  is  chosen  such  that  the  spectral  plane  wave  either  propagates  or  decays 

exponentially  as  it  moves  away  from  the  xy  plane.  In  general,  for  a  lossless  host 

medium  (Im(fco)  =  0)  there  will  be  a  finite  number  of  propagating  waves  and  an 

infinite  number  of  decaying  evanescent  waves. 

The  incremental  magnetic  field  is  found  simply  as 


1  Idl'  e-j*o*  *± 

dH  =  ±vxdA  =  -^-TEE - 

no  2duvudv  ~  ~  r. 


(2.52) 


and  the  incremental  electric  field  is  found  as 

1  Idl'i jo  ^  c-itoRf± 

dE  = - V  x  dH  = - —  >  V - e+ 

iu,eoV  2dvVyd.it  %  r,  ** 


*  where, 


(2.53) 


e±  =  (a  x  ?±)  x  f±  =  (f±  •  a)r±  -  a. 


(2.54) 


current 
element 
segment  nj 


nj 


Figure  2.5:  Vector  relationship  for  the  center  Q  =  0  current  segment  nj. 


If  the  reference  point  moves  from  the  origin  to  R',  the  resulting  field  is  the  same 
as  if  the  reference  point  stays  at  the  origin,  and  the  field  point  moves  to  R  —  R'.  In 
this  manner,  dE  becomes 


The  total  electric  field  of  a  linear  current  element  nj,  denoted  as  Fnj(V),  is  now 
determined  by  integrating  Equation  (2.55)  along  the  length  of  the  current  element. 
Typically,  a  current  element  nj  will  be  monopole  j  (j  =  1  or  2)  of  expansion  function 
n.  Note  that  in  Equation  (2.55)  the  only  variation  on  the  source  point  R'  is  in  the 
last  exponential  term.  As  shown  in  Figure  2.5,  if  current  element  nj  is  straight,  the 
source  point  R'  can  be  written  as 


R^R^  +  fin,!'  for  0<l'  <dnj  (2.56) 

where  Rnj  is  the  position  vector  from  the  origin  to  the  reference  end  of  current  element 
nj,  and  d„j  is  the  length  of  segment  nj. 

Figure  2.6  shows  that  there  are  three  possible  distinct  regions  that  the  field  point 
R  can  be  within.  The  field  in  Region  (I)  consists  only  of  left-going  waves,  so  the  f_ 
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spectral  vector  is  chosen  and  the  integration  is  over  the  entire  range  of  the  current 
element.  Similarly,  the  field  in  Region  (III)  consists  only  of  right-going  waves,  so  ?+ 
is  chosen  and  the  integration  is  over  the  entire  current  element.  In  this  manner,  the 
electric  field  in  Region  (I)  can  be  written  as 


2duvvdv 


EE 


®nj-  Pnj-  (0»  hi  ) 


(2.57) 


and  the  electric  field  in  Region  (III)  can  be  written  as 


enj+  dnj)' 


(2.58) 


The  function  ,  b)  is  termed  the  njtk  source  pattern  factor  and  is  defined  as 

P„‘J±(a,4)  =  £  (2.59) 

where  Fnj(V)  is  the  current  distribution  along  the  length  of  current  segment  nj. 
The  source  pattern  factor  is  evaluated  in  greater  detail  in  Appendix  A.  In  Region 
(II)  the  electric  field  at  point  R  consists  of  left-going  waves  from  that  section  of 
the  current  element  denoted  by  nj-  and  right-going  waves  from  that  section  of  the 
current  element  denoted  by  nj+.  Thus,  the  field  at  point  R  in  Region  (II)  is  the  sum 
obtained  by  choosing  the  r_  spectral  vector  and  integrating  Equation  (2.55)  over 
nj-  plus  choosing  the  r+  spectral  vector  and  integrating  Equation  (2.55)  over  nj+. 
Therefore,  the  electric  field  in  Region  (U)  can  be  written  as 

e£(r)  =  mtzEE  — ; - (<■»,+ a,+)  + 

ZOuVyOy  nv  n>  r. 


H  “  en>-  Pnj-  (®"i- » bnj- )  •  (2.60) 

The  limits  of  integration  are  from  a„j+  to  bnj+  for  evaluating  the  P*J+  integral  and 
from  Onj-  to  bnj-  for  evaluating  the  P’j_  integral.  These  limits  depend  upon  the 
specific  wire  segment  geometry,  and  are  given  in  Section  A.l. 

Note  that  if  the  current  element  is  in  a  plane  parallel  to  the  xy  plane,  then  Region 
(II)  vanishes  and  there  are  only  two  regions.  In  this  case,  Equations  (2.57)  and  (2.58) 
can  be  used  to  determine  the  electric  field  everywhere. 
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(2.61) 


In  a  similar  manner,  the  magnetic  fields  in  the  three  distinct  regions  are: 

H®(R)  =  rr—j-EE - r - (**  x 

la^Vydy  nu  n>  r, 

/Till  1  . _ .  __ 

Hi,n(R)  =  5T^-rEE£ - r - (*.,  X  (2.62) 

zOuVyOv  B>  n<  r, 

H^’fR)  =  j^EE  [ - - - («.,  x  f+)P„V+KJ+.6»J+)  + 

+ - - - (**  X  §_)/J_(«*_,fc*_)  .  (2.63) 

It  should  be  noted  that  the  fields  evaluated  in  this  section  are  for  a  single  current 
element  nj.  This  is  adequate  for  evaluating  the  fields  from  a  monopole  expansion 
function.  However,  a  dipole  expansion  function  is  composed  of  two  current  segments, 
and  the  procedure  of  this  section  must  be  repeated  for  each  section.  As  shown  in 
Figure  2.5  and  in  Equation  (2.56),  the  vector  Rnj  always  points  from  the  origin  to 
the  zero  current  endpoint  of  current  segment  nj,  and  the  unit  vector  a*,-  points  from 
the  zero  current  endpoint  to  the  unity  current  endpoint. 

2.5  Fields  of  a  3D  Volume  Array  of  Current  Ele¬ 
ments 

This  section  uses  the  results  of  Section  2.4  to  evaluate  the  exact  electric  and  magnetic 
fields  of  a  3D  volume  array  of  linear  current  elements  [32,  Gh.  2].  As  shown  in 
Figure  2.7,  the  3D  volume  array  consists  of  an  infinite  parallel  stacking  of  the  2D 
planar  array  of  elements  in  the  host  medium  (po,  eo).  The  2D  layers  are  referenced  by 
the  index  for  —  oo  <  K,  <  oo,  and  are  spaced  a  distance  of  <£*,  in  the  unit  direction 
w  =  Wgk  +  wvy  +  wxi.  Therefore,  the  vector  from  the  reference  point  on  the  fc*,  =  0 
plane  (the  origin)  to  the  equivalent  point  on  any  fc*,  ^  0  plane  can  be  skew  to  the  xy 
plane,  and  is  given  as  w.  The  current  elements  on  each  plane  will  be  weighted 
by  the  plane  wave  propagating  through  the  volume  array.  The  current  weighting 
corresponding  to  the  plane  wave’s  variation  in  the  u  and  v  directions  have  already 
been  included  in  the  analysis  of  Section  2.4.  The  current  weighting  corresponding  to 


Figure  2.7:  A  3D  volume  array  of  made  of  stacked  2D  planar  arrays, 
the  plane  wave’s  variation  in  the  w  direction  is 

C{kw)  =  =  e->(fc«d*)  where 

dig  =  d|o(kex  Wx  "f  key  ^y  d"  ke*  ®i). 

Primarily  of  interest  are  the  fields  in  the  center  lattice  cell  about  the  origin.  These 
fields  are  the  sum  of  all  the  2D  planar  element  fields  for  all  planes  —  oo  <  K>  <  oo. 
The  fields  of  the  k*,  =  0  planar  array  were  determined  in  Section  2.4.  The  fields  of 
the  planes  referenced  by  k*  =  +1,  +2,  +3,  •  •  • ,  +oo  consist  of  only  left -going  waves 
in  the  center  lattice  cell.  Thus,  E^(R)  and  H$(R)  of  Equations  (2.57)  and  (2.61) 
must  be  weighted  by  C(kw)  and  summed  over  positive  k».  Similarly,  the  fields  of 
the  k«,  =  —1,  —2,  —3,  •  •  • ,  —  oo  planes  consist  of  only  right-going  waves  in  the  center 
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lattice  cell.  Thus,  Ej^'(R)  and  H^^(R)  of  Equations  (2.58)  and  (2.62)  must  be 
weighted  by  C(kw)  and  gummed  over  negative  k^.  It  should  be  noted  here  that  the 
z  range  spanned  by  elements  in  adjacent  At*  planes  do  not  overlap.  (See  Figure  2.6.) 
This  requirement  insures  that  there  are  no  Region  (II)  type  fields  from  planes  where 
For  both  positive  and  negative  to  account  for  the  shift  in  position  to  the 
plane  indexed  by  At*,,  the  vector  R nj-  (of  the  field  equations  for  Regions  (I)  and  (III)) 
must  be  replaced  by  +  k^d*, w. 

In  this  manner,  the  total  electric  field  in  the  center  lattice  cell,  due  to  the  3D 
volume  array  of  all  the  njih  current  elements,  can  be  written  as 

E5OT(R)  =  E^(*„  =  0,a)  + 

+  £  <7(M E#(*.,R)  +  2  C(fc)E<J“>(*„,R)-  (2.64) 


*»=+! 


*.=-1 


The  electric  field  E;A-  (fc*,,R)  can  be  written  as 


enj—  Pni-  (Oi  )•  (2.65) 


Inserting  E^(^,  R)  into  the  first  summation  of  Equation  (2.64)  and  rearranging  the 
summation  order,  it  is  obtained  that 

2  C(*»)E$(*»,R) 


*,=+1 


Vo 

2{2||Vyl 


E  5 - ; - «*-  P„V(0.<,).  (2.66) 

*„=+!  /  T* 


Similarly,  the  second  summation  of  Equation  (2.64)  becomes 

2  0(4.)  E<Jn>(*.,R) 


*»  =  -! 


Vo 

2dnVy 


E  «•**>)*  r  (2.67) 

*»=-!  /  1 


For  Equations  (2.66)  and  (2.67),  0±  =  d*  —  fcod*,(w  •  r±)  =  d«,[ke  —  fcor±]  •  w. 
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The  innermost  summation  on  h*  can  be  evaluated  in  closed  form  by  making  use 
of  the  summation  identity 


tel 


1  -  e-tf 

where  (3  is  a  complex  number.  The  summation  on  h*  in  Equation  (2.66)  is  a  straight¬ 
forward  application  of  this  identity.  However,  the  summation  on  h*  in  Equation  (2.67) 
requires  a  change  in  the  summation  index  of  k*  — »  —  k'wt  such  that  the  summation  is 
over  positive  k'w.  The  results  are 


+oo 

E 

*«=+! 


_  e_ 

1  -  e~i&~ 


and 


S»=-1  1 


e#+ 


(2.68) 


Combining  these  results,  the  total  electric  field  in  the  center  lattice  cell  can  be  written 


as 


E™T(R)  =  E„)(fc.=0,R)  + 

/ 

fin  e-rt-  \ 

+  [(l -«-*-)  7,  e">-^-(°'*»)+ 

__j  -  eni+PV+(0,dnj)  . 


(2.69) 


Performing  a  very  similar  operation  on  the  magnetic  field,  the  total  magnetic  field  in 
the  center  lattice  cell  can  be  written  as 


Hj0T(R)  =  HnJ(fc*  =  0,  R)  + 


+ 


2duVydv 


- - - -(•„,  X  f-)P'i_(0,4i)  + 


+  ^  _  ej0+)  (•*«'  X  *+)^*»+(®»^i)  • 

For  simplicity,  is  is  worthwhile  to  write  these  total  fields  as 

e£0T(r)  =  =  0,R)  +  Eni(lt«,  >  0,R)  +  EBi(Jb«  <  0,R) 

hT°t(R)  =  Hni(fc„  =  0,R)  +  Hni(kw  >  0,R)  +  IMfc*  <  0,R). 


(2.70) 


(2.71) 

(2.72) 
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Chapter  3 

Evaluation  of  the  Impedance 
Matrix 


Typical  elements  of  the  impedance  matrix  are  given  by  Equation  (2.10)  for  the  Zmn 
terms,  and  by  Equation  (2.11)  for  the  A Zmn  terms.  The  evaluation  of  the  A Zmn 
terms  is  straightforward  and'  is  presented  in  Section  3.1.  The  Zmn  terms  are  more 
complicated  to  evaluate,  as  will  be  outlined  next. 

Recall  that  expansion  function  n  can  be  either  a  dipole  (j  =  1  and  2)  or  monopole 
(j  sss  1)  expansion  function.  Similarly,  weighting  function  m  can  be  either  a  dipole  (» 
=  1  and  2)  or  monopole  (*  =  1)  weighting  function.  Therefore,  there  are  four  mutual 
impedance  possibilities.  Following  the  polarity  convention  presented  in  Section  2.2.1, 
the  four  mutual  impedance  possibilities  are: 


1.  *  =  2  ,j  =  2: 

2.  i  =  l,j  =  2: 

3.  i  =  2,  j  =  1: 

4.  »  =  1,  j  =  1: 


2*n  —  Zmi<ni  Zml,n2  ^r»2,nl  +  Zm  J,n2, 

Zmn  —  Zml,nl 

Zmn  =  Zml,nl  ~  ^m2,nl>  and 
Zmn  ~  Zm  1,1,1 


where  Zmi>nj  is  the  mutual  impedance  between  monopole  j  of  expansion  function  n 
and  monopole  s  of  weighting  function  m. 

From  Equation  (2.10),  on  a  monopole  to  monopole  basis,  it  is  seen  that 


=  E  cQ4?< 

Q 


mtjij 


(3.1) 
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where 


zS«  =  -/  (3.2) 

JVmi 

in  which  E|J  is  the  electric  field  of  jj^-,  and  the  integration  is  over  K,,,  the  volume  of 
monopole  *  of  weighting  function  m.  The  electric  field  EjJ  results  from  both  the  axial 
and  the  radial  components  of  given  in  Equations  (2.15)  and  (2.16).  Therefore,  it 
is  convenient  to  write  E^  as 

Eg  =  Eg'  +  Kg’  (3.3) 

where  E^1  is  the  electric  field  of  the  axial  current  component  and  is  the  electric 
field  of  the  radial  current  component.  Similarly,  it  is  convenient  to  write  Z^inj  u 


_  pQs  .  «rQ p 


Z&j  =  -i  Eg'  .  Wml  dv 

.j  *  “mi 

Z**?  -  -  f  E*^  •  W„,  dv. 

*Mni 

The  evaluation  of  Z^fnj  is  presented  in  Section  3.2. 

The  Z%?nj  terms  are  evaluated  for  all  values  of  Q  at  once  by  noting  that 


Q 


EC*# 

"»•  I  Q 


Wmi  dv 


where 


Ecqe£ 

Q 


=  bT0T 


n J 


(3.4) 

(3-3) 

(3-6) 


(3.7) 


with  given  in  Equation  (2.69)  or  (2.71).  The  evaluation  of  Equation  (3.7)  is 

rather  involved,  and  is  presented  in  Section  3.3. 


3.1  Evaluation  of  A Z 

A  typical  A Zmn  term  is  given  by  Equation  (2.11),  which  involves  only  the  expansion 
and  weighting  functions  defined  in  the  center  Q  =  0  lattice  cell.  Note  that  the 
A  Zmn  term  will  be  non-sero  only  when  expansion  function  n  overlaps  with  weighting 
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monopole  i 

weighting  function  m 

monopole  j 
expansion  function  n 


monopole  t 

weighting  function  m 
monopole  j 

expansion  function  n 


Figure  3.1:  The  two  possible  monopole  to  monopole  overlap  arrangements. 

7 

function  m.  Since  the  expansion  functions  are  made  up  of  monopoles,  the  overlapping 
is  treated  on  a  monopole  to  monopole  bans.  Figure  3.1  shows  the  two  possible  ways 
that  monopole  j  of  expansion  function  n,  can  overlap  with  monopole  t  of  weighting 
function  m,  for  i,j  =  1  and/or  2. 

The  A  Z  contribution  when  monopole  j  of  expansion  function  n  overlaps  with 
monopole  i  of  weighting  function  m  can  be  written  as 

r  /; 

where  d  is  the  length  of  the  overlapping  monopole  segment.  Noting  the  expressions 
for  the  expansion  and  weighting  functions  in  Section  2.2.1,  Equation  (3.8)  becomes 

=  M  -  eo)*""™  (3-9) 

where  Fm^„j  depends  upon  the  axial  variations  of  the  overlapping  segments,  and  is 


written  as 


F m»,nj  —  j  F *)  F nj(ho  *)  dx. 
JO 


(3.10) 
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In  both  overlapping  cases  of  Figure  3.1,  the  axial  variation  for  the  weighting 
monopole  mi  is 

d 


w..).f*W*,s*s' 

0  otherwise. 


(3.11) 


If  the  overlapping  geometry  is  as  shown  in  Figure  3.1(a),  then  the  axial  variation 
of  the  expansion  monopole  nj  is  the  same  as  for  the  weighting  monopole  mi,  and 
Equation  (3.10)  becomes 


_  hp  d  —  cos(feo  d)  sin(fco  d) 
mi  ni  2  ho  «mJ(ho  d)  ' 


(3.12) 


If  the  overlapping  geometry  is  as  shown  in  Figure  3.1(b),  then  the  axial  variation  of 
the  expansion  monopole  nj  is 


if0<z< 
0  otherwise, 


(3.13) 


and  Equation  (3.10)  becomes 


„  sin(hod)  —  hod  cos(hod) 

— 2*oBn4(*»  i)  ’ 


(3.14) 


The  minus  sign  in  Equations  (3.13)  and  (3.14)  account  for  the  fact  that  the  current 
on  the  monopoles  in  Figure  3.1(b)  are  of  opposite  polarity. 


3.2  Evaluation  of  Z**p 

This  section  evaluates  given  in  Equation  (3.6).  As  indicated  in  Section  2.2.1, 
the  MM  expansion  functions  have  a  radial  current  component,  defined  on  a  per 
monopole  basis  and  given  by  Equation  (2.16).  However,  the  MM  weighting  functions, 
Wm„  given  by  Equations  (2.20)  and  (2.21),  contain  only  an  axial  current  component. 
The  pnj-directed  current  J£j  will  produce  an  electric  field,  which  will  be  highly  lo¬ 
calised  to  an  axial  field  along  its  centerline.  Therefore,  Z^nj  will  be  approximated 
as  sero,  unless  if  monopole  j  of  expansion  function  n  overlaps  with  monopole  *  of 
weighting  function  m.  In  this  manner,  only  the  radial  currents  in  the  center  Q  =  0 
wire  object  contribute  to  the  impedance  matrix.  In  summary,  =  0  for  all  Q  /  0 
and  for  all  nonoverlapping  segments. 
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As  illustrated  in  Figure  3.1,  and  explained  in  Section  3.1,  the  overlapping  is  treated 
on  a  monopole  to  monopole  basis.  If  monopole  j  of  expansion  function  n  overlaps 
with  monopole  i  of  weighting  function  m,  then  the  contribution  to  the  impedance 
matrix  element  resulting  from  Jnjy  is  written  as 


=  jf  ff  E£  y/~pdpd*d2  (3.15) 


where  d  is  the  length  of  the  overlapping  monopole  segment,  and  is  the  electric 
field  of  J*j.  Referencing  the  work  of  Richmond  and  Newman  [33],  this  becomes 


7* 


C[l-Jo(fc,q)] 
jwe ohg  mini 


(3.16) 


where 

=  jf  Ku *)  *)  i*.  (3.17) 

The  prime  '  denotes  differentiation  with  respect  to  z,  and  F^koz)  and  Fnj(koz)  are 
given  by  Equations  (3.11)  and/or  (3.13),  depending  on  the  overlap  geometry.  If  the 
overlap  geometry  is  as  in  Figure  3.1(a),  then 

i  <0  co.(*„ i)\  (3.18) 

and  if  the  overlap  geometry  is  as  in  Figure  3.1(b),  then 


Ki,nj  =  2siJlk0d)  ^  d  COS^*°  ^  +  “n(*°  ^ '  (3  19) 

The  opposite  polarity  of  the  overlap  case  of  Figure  3.1(b)  has  been  included  into 
Equation  (3.19). 


3.3  Evaluation  of  Z* 


This  section  evaluates  the  monopole  to  monopole  impedance  given  in  Equa¬ 
tion  (3.7).  In  general,  the  scalar  product  of  with  Wmt-  is  integrated  throughout 

Vmi,  the  volume  of  monopole  *  of  weighting  function  m.  Thus,  the  impedances  require 
a  triple  volume  integration  throughout  However,  in  all  cases  this  can  be  reduced 
to  a  single  integration  in  the  axial  direction  of  the  weighting  function  monopole  mi. 
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If  the  monopoles  mt  and  nj  do  not  physically  touch  or  overlap,  then  the  currents 
can  be  approximated  as  line  sources  of  the  same  axial  current  variation  and  strength, 
located  at  the  centerlines  of  the  monopole  segments.  In  this  case  no  stagger  of  the 
weighting  monopole  is  necessary.  However,  if  both  monopoles  are  in  the  same  plane, 
and  that  plane  is  parallel  to  the  xy  plane,  then  a  slight  stagger  of  the  weighting 
monopole  out  of  this  plane  can  accelerate  the  convergence  of  the  spectral  summations 
involved  in  the  impedance  evaluation.  This  is  due  to  the  slight  exponential  decay  of 
the  evanescent  waves  in  the  spectral  summation  of  plane  waves. 

If  the  monopoles  mi  and  nj  do  physically  touch  or  overlap,  then  the  triple  volume 
integration  is  reduced  to  a  single  integration  by  employing  the  method  of  equivalent 
wire  radius,  as  used  by  Newman  for  thin  material  wires  [25].  Basically,  this  results  in 
the  elimination  of  the  dp  and  dj>  integrations  by  staggering  the  weighting  monopole 
an  equivalent  wire  radius  off  its/  centerline.  The  theory  of  this  equivalent  wire  radius 
is  presented  in  Appendix  B,  and  the  rules  for  staggering  the  weighting  monopole  mi 
are  given  in  Appendix  C. 

When  integrating  in  the  axial  direction  of  the  monopole  mt,  it  is  worthwhile  to 
note  that  the  position  vector  to  a  point  on  monopole  mt  can  be  written  as 

R  =  Rn,,  +  a mil  iorO<l<dmi  (3.20) 

where  R*,,  is  the  position  vector  from  the  origin  to  the  reference  end  of  monopole  mt, 
and  dmi  is  the  length  of  monopole  mt.  Recall  that  for  the  PWS  monopoles  used  in  this 
paper,  the  reference  end  is  the  zero-current  end.  Equation  (3.20)  is  a  direct  result 
of  Figure  2.5  applied  to  monopole  mt  instead  of  monopole  nj,  where  the  position 
variable  l'  has  been  replaced  by  l. 

The  impedance  Z^i>ni  involves  the  electric  field  of  the  3D  array  of  current 

segments  nj.  is  given  in  Equation  (2.71),  where  it  is  represented  as  contribu¬ 

tions  from  the  three  ranges  of  the  lattice  summation  index  h*,  (fc»  =  0 >  0  and 
k*,  <  0).  Thus,  it  is  convenient  to  write  Z^in-  as 

ZA ,»  =  ZA*„,  +  Z£,y  +  ZAU  ’  (3-21) 

-  -  jf  EW(*.=0,R)-Wm(*  (3.22) 

*'n»i 
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(3.23) 

(3.24) 


ZZnj  =  -  /  >  o,  R)  •  Vfmi  dv 

KU,  =  -  }Vmt  E.,(k.<0,R)-W.,J.. 

In  all  the  above  cases,  the  weighting  monopole  becomes  the  filament  source 

Wm,  =  fi™  (3.25) 

3.3.1  The  Evaluation  of 

The  form  of  =  0,  R)  used  in  the  evaluation  of  Z^i  nj  depends  upon  the  geo¬ 

metrical  arrangement  of  monopoles  mi  and  nj.  Figure  D.l  and  Appendix  D  describe 
the  eight  possible  geometrical  arrangements,  which  are  labeled  as  Case  1  -  Case  8. 
Basically,  Equations  (2.57), (2.58)  and/or  (2.60)  are  used  to  evaluate  Enj(fc*  =  0,R) 
in  Regions  (I),  (II)  and/or  (III),  as  required  by  the  geometry  of  monopoles  mi  and 
nj.  The  eight  possible  cases  for  the  evaluation  of  Z£-nj  are  analyzed  here. 

Case  1  and  Case  2: 

Cases  1  and  2  are  the  simplest  because  the  electric  field  of  expansion  monopole  nj 
across  weighting  monopole  mi  consists  solely  of  left-going  or  right-going  waves,  re¬ 
spectively.  For  Case  1,  Enj(A»  =  0,R)  =  E$(R)  of  Equation  (2.57),  and  for  Case 
2,  E nj(kw  =  0,R)  =  E^*(R)  of  Equation  (2.58).  Since  monopole  mi  is  entirely 
within  either  Region  (I)  or  Region  (III)  of  expansion  monopole  nj,  only  one  expres¬ 
sion  for  Enjffcu  =  0,R)  is  required,  and  the  limits  of  integration  are  over  the  entire 
length  of  monopole  mi  in  that  region.  The  only  spatial  variation  of  E nj(K>  —  0,  R)  is 
the  exponential  term  e~jtoR  t±  inside  the  double  spectral  summation  of  plane  waves. 
Thus,  the  integration  only  involves  this  exponential  term  and  the  weighting  function 
variation  of  Equation  (3.25). 

Inserting  Enj(fc«,  =  0,R)  and  Wm,  into  Equation  (3.22)  and  integrating,  it  is 
obtained  for  Case  1  that 

X  (e*,.  .  imi)  Ptj-i 0,  dn;)P;,-(0,  <L,)  (3.26) 
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and  for  Cute  2 


-  2duvvd„  £ 


EE - ; - x 


X  (e»i+  '  *"M )  ^nj+  (® » )^m»+  (®» ^I»t )  (3.27) 

where  P,^l±(a,6)  is  termed  the  m***  receive  pattern  factor  and  is  defined  as 

P;.±(«,4)  =  £  F..(0  «-'*•'(*-  '*)  dl  =  Cl*  [AU(6)  -  4.1±(.)]  (3.28) 

where 

CL±  =  MH(*0  [1  ~  (*-.  ■  f*)1]  (3'29) 

A^,±(*)  =  •  ?±)  «in(*o  *)  -  cos(*o  *)]  •  (3.30) 


Case  3: 


Case  3  is  more  complicated  because  Bnj(kw  =  0,R)  =  E®  (R)  of  Equation  (2.60), 
which  consists  of  both  right-going  and  left-going  waves  across  weighting  monopole 
mi.  However,  monopole  mi  is  entirely  within  Region  (II)  of  expansion  monopole 
nj,  so  the  limits  of  integration  are  over  the  entire  length  of  monopole  mi  using  the 
Region  (II)  field.  The  spectral  summations  contain  the  usual  spatial  variation  of  the 
exponential  term  e-i*o  Moreover,  the  source  pattern  factors  Pnj±(*nj±iKj±) 


contain  a  spatial  variation  due  to  the  z  dependence  of  their  integration  limits,  given 
in  Section  A.l.  With  these  considerations  for  Case  3,  Z^^-  becomes 


71=  _ 


Vo 


2duVvdv  X 


EE 


r  g— [H-mi — 


n. 


(0,  dmi  )+ 


"f"  "  Gtni.nj—  (0,  dmi ) 

Tt 

where 


(3.31) 


Gm,,„j±(a,  b)  =  £ (e.i±  •  (3.32) 

Note  that  Pnj±(anj±,bnj±)  has  an  l  dependence  due  to  the  source  integration  limits. 
The  scalar  product  (enj+  -a™,)  appears  inside  the  integration  so  that  open-current  end 
charge  contributions  can  be  isolated  and/or  removed  (see  Section  A.2).  Gmi<nj±(a,b) 
is  a  rather  involved  integration  and  is  evaluated  in  Appendix  E. 
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Case  4,  Case  5  and  Case  0: 


Note  from  Figure  D.l  that  in  Case  4  the  weighting  monopole  is  partially  in  Region  (I) 
and  partially  in  Region  (II)  of  expansion  monopole  nj.  Therefore,  E =  0,R)  = 
E$(R)  of  Equation  (2.57)  when  integrating  over  the  portion  of  monopole  mi  that 
is  in  Region  (I),  and  Enj(fc»  =  0,R)  =  E®(R)  of  Equation  (2.60)  when  integrating 
over  the  portion  of  monopole  mi  that  is  in  Region  (II). 

Similarly,  for  Case  5  the  weighting  monopole  is  partially  in  Region  (III)  and 
partially  in  Region  (II)  of  expansion  monopole  nj.  In  this  case,  E„s(kw  =  0,R)  = 
EjJ^(R)  of  Equation  (2.58)  when  integrating  over  the  portion  of  monopole  mi  that 
is  in  Region  (III),  and  E„j(kw  =  0,R)  =  E^(R)  of  Equation  (2.60)  when  integrating 
over  the  portion  of  monopole  mt  that  is  in  Region  (II). 

When  integrating  in  Region  (I)  or  Region  (HI),  the  spatial  variation  of  Enj(fc„,  = 
0,R)  consists  only  of  the  exponential  term  However,  when  integrating  in 

Region  (II),  the  spatial  variation  of  E„j(kw  =  0,  R)  includes  the  exponential  term 
and  the  source  pattern  factor.  With  these  considerations,  for  Case  4  Z^~  nj  becomes 


z*=  .  =  — 5° — yy 

2duvvdv  ~  “ 

~  (e»i-  *  ®mi)  Pnj-  dnj)Pmi-{a™iubmii)  + 

"i 

Gfnt.nj—  (®mi2)  "(■ 

T* 

g j ^0  [®»m«  —  1 

H  G m« ,nj+  ( Omi! »  ^mt'2 ) 

r* 

and  for  Case  5  Z£-nj  becomes 


(3.33) 


7*=  =  _  *7°  V  V 

w<>ni  2 

g “  jh)  [B-m » “  j  ] + 

~  (ei»i+  *  ®m»)  (°mi3»  &m«3)  + 

Tz 

“1"  Gmi>nj—  (fltrn'2)  ^mij)  “f* 

** 

H - Gmi  ,n;+  (®mi2  j  ^mi'2)  I  • 


(3.34) 
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In  either  case,  the  limits  of  integration  (am,i ,  2>m,i ),  (cw,  )  and  (0^3,  Jwj)  repre¬ 
sent  the  portion  of  weighting  monopole  mi  that  is  in  Region  (I), (II)  or  (III),  respec¬ 
tively.  These  limits  are  defined  in  Appendix  F. 

From  Figure  D.l,  it  is  seen  that  Case  6  is  an  extension  of  either  Case  4  or  Case 
5  in  that  monopole  mi  extends  over  Regions  (I), (II)  and  (III)  of  monopole  nj.  The 
same  statements  made  above  for  Cases  4  and  5  hold  true  for  Case  6,  so  that  Z£~nj  is 
written  as 


7*=  _ 


Vo 

2^utJy(£j/ 


EE 

fly  fly 


(enj-  *  •mi )  ^nj - (®»  )^m»- (®mil >  ^rail )  + 

T* 

g““  j^O  "" 

H  (®nj+  *  •mi)  ^>nj'+(0)  ^nj)-Pmi-(ami3j  bmrt)  + 

*** 

g ~ j&O  [H*mi  ” j  ] 

*f  1  Grni%nj-  (<*mi2f  4” 

g"A  [R«m»  “®t»j  ]*^ + 

4 - Gmi  ,n}+(^ffli2i^nu2)|  • 

The  limits  of  integration  are  given  in  Appendix  F. 


(3.35) 


Case  7  and  Case  8: 


From  Figure  D.l,  it  is  seen  that  Case  7  can  be  viewed  as  a  special  case  geometry  of 
Case  6,  where  Region  (II)  of  monopole  nj  has  vanished.  Thus,  Z%~nj  consists  of  only 
the  Region  (I)  and  Region  (III)  contributions,  and  is  written  as 


Z*=  = 


Vo 


2dvvudv 


EE 

n*  n« 


rx 


(®nj-  *  •mi)  l*nj-(®)^»j)^ mi-(amilj  ^mil)  + 


+ 


g"‘i^oI^mi"^»»i],^+ 


Tt 


(enj+  •  •mi)  ■Pn,>(0»d»»i)-Pm<+(Omi3»6m,-3) 


(3.36) 


Case  8  is  similar  to  Case  3,  except  that  in  Case  8  monopole  mi  does  not  have  any 
extent  in  the  z  direction.  However,  it  is  still  entirely  within  Region  (II)  of  monopole 
nj.  The  limits  of  integration  become  independent  of  z,  and  the  integrations  of  (3.32) 
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no  longer  are  nested.  For  case  8,  611X1  be  given  by  Equation  (3.31),  except 

that  the  pattern  factor  integral  term  now  simplifies  to 


Gmi,txj±  (®»  6)  —  Pmi±  (al  Pnj±  (°»i±  »  ^*j±  ) 


(3.37) 


where  the  limits  of  integration  are  given  in  Appendix  F. 

3.3.2  The  Evaluation  of  and  Z£nj 

The  electric  field  from  all  planes  indexed  by  fc*,  >  0  consists  only  of  left-going  plane 
waves  in  the  center  Q  =  0  lattice  cell,  and  is  expressed  as 

fin  /  \  I®-— •*— 

Bn'(*"  >  °’R)  =  2d^??(l~e-^-j  7, 

(3.38) 

The  only  spatial  variation  of  Enj(fc*,  >  0,R)  is  the  exponential  term  e~ik°R  t-  inside 
the  double  spectral  summation  of  plane  waves.  Employing  Equation  (3.20),  Z£$nj- 
can  be  expressed  as 


ZI*  = 


Vo 


^  ^  f  e->0-  ^  j-iMlni-lajl't- 


mi,nj  ~  2duVydv 

X  (enj-  *  dmi)  ^j_(0,  dnj)Pm|_(0,dml) 

where  •P,Jri_(0,dmi)  is  given  in  Equation  (3.28). 

Similarly,  Z^i<nj  can  be  expressed  as 

ijq _ (  er&*  \ 


(3.39) 


Z'<  .  =  • 
w  2duvydv 


X  (enj+  '  ■tni)^>Bi+(®»^"j)^mi+(®><^»»*) 


(3.40) 
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Chapter  4 

Average  Eigenfunction  Electric 
and  Magnetic  Fields 


This  chapter  presents  the  evaluation  of  the  electric  and  magnetic  eigenfunction  fields 
averaged  over  the  center  Q  =  0  lattice  cell.  These  fields,  denoted  as  (Ea,H#),  are 
defined  as  ' 

E‘=^X.E(R>A'  (“i 

H* = i;  L  (4-2) 

where  Av  is  the  volume  of  a  lattice  cell  in  the  3D  array  and  V°  represents  the 
integration  limits  for  the  center  Q  =  0  lattice  cell.  For  the  evaluation  presented  here, 
the  3D  lattice  array  consists  of  perpendicular  axes,  i.e., 

U  =  X,  V  =  y,  W  =  £, 

du  —  dxt  dy  —  dy ,  dy,  =  dt 

and  therefore  Av  =  dx  dy  dx. 

The  eigenfunction  fields  (E,  H)  are  the  fields  radiated  by  the  eigenfunction  cur¬ 
rents,  so  for  the  PMM  solution 

N  1(2) 

E(R)«  (4-3) 

n=l  j=l 
N  1(2) 

H(R)  «  E  J.  E  (-1)  Hkj(R)  iv  (4.4) 

nal  j=l 

where  the  summation  on  j  is  j  =  1  for  a  monopole  and  j  =  1,2  for  a  dipole  expan¬ 
sion  function.  The  sign  factor  (— 1)*'+1  accounts  for  the  polarity  convention  of  the 
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monopoles  making  up  the  expansion  functions.  In  evaluating  the  fields  of  a  given  MM 
expansion  function,  the  fields  of  the  radial  current  component  are  neglected,  and  thus 
Equations  (2.71)  and  (2.72)  are  used  to  evaluate  the  fields  of  current  segment  njf. 
Noting  that  Equations  (2.71)  and  (2.72)  represent  the  fields  as  contributions  from 
three  ranges  of  the  lattice  summation  index  (k«,  =  0,fc»  >  0  and  <  0),  the 
averaged  eigenfunction  fields  are  written  as 
i  n  i(a) 

E*  *  £  /.  £ (-1) <+'  =  0)  +  $(*.  >  0)  +  I*  (i.  <  0)]  (4.5) 

n=l  j=l 
1  N  1(3) 

H°  »  T-  £  A.  E  =  «)  +  >  0)  +  <  »)]  (4.6) 

n=l  j=l 

where 

$(*.  =  o)  =  jT ,  E„(4»  =  0,R)A. 

<5<*-  > 0)  =  /v,  E.y(k»  >  0.R)*. 

<  0)  =  jvt  E„,(fc.  <  0,  R)  dv 
l3(*.  =  «)  =  X,  H.y(fci.  =  0,R)«fo 
>  0)  =  X,  Hn)(*»  >  0, R)  (it; 

Ig(*.  <  0)  =  X.  H”i(*-  <  ®. R)  *•  (4-7) 

Recall  from  Figure  2.6  that  in  general  there  are  three  possible  regions  for  a  field 
point  in  the  center  Q  =  0  lattice  cell,  due  to  the  K,  =  0  array  of  current  segments 
nj.  These  regions  are  defined  by  the  z  values  of  the  current  segment  endpoints. 
Figure  4.1  shows  the  limits  on  the  z  integration  corresponding  to  the  three  regions. 
It  is  convenient  to  write 

£,(*.  =  o)  =  ^  +  £!n)  +  #ra) 

O*- = °) = C(I)  +  i5(n>  +  J5(m) 

where 

£fI) = X i  ■$<*>* 


(4.8) 

(4.9) 
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width  of  lattice  cell 


Figure  4.1:  The  z  integration  limits  for  evaluating  the  average  eigenfunction  fields. 

£(UI)  =  X„I 

*3(I) = L 

>3(n)  =  Xn 

C(III)  =  Xra  H*V>‘h’-  («•">) 

where  the  integration  limits  V^,  and  refer  to  the  Region  (I), (II)  and  (III) 
portions  of  the  center  lattice  cell  for  the  Ki  —  Q  current  segment  nj,  respectively.  The 
electric  fields  E$(R),  E^R)  and  E^II}(R)  are  given  by  Equations  (2.57), (2.60) 
and  (2.58),  respectively,  and  the  magnetic  fields  H$(R),H^(R)  and  H^^(R)  are 
given  by  Equations  (2.61), (2.63)  and  (2.62),  respectively.  Note  that  if  current  segment 
nj  is  in  a  plane  parallel  to  the  xy  plane,  then  *3  =  *1,  and  Region  (II)  vanishes,  and 
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4.1  kw  =  0,  Region  (I) 


From  the  analysis  in  Sections  2.4,  it  was  found  that  for  Region  (I)  and  Region  (III), 
the  only  spatial  variation  of  the  fields  is  the  exponential  term  e~J*°*'**  inside  the 
double  spectral  summation  of  plane  waves.  Thus,  ijy  ^  and  1^^  can  be  expressed 

C(I) - <«»> 

1°(1)  =  X  (4.12) 

^“x“v  n«  n.  r*  JVl 

The  limits  of  integration  for  Region  (I)  are 

\  <  *  <  *t 


V1  =  J  _i  < 


s  y  s 


<  i 

—  3 


-7  <  *  < 


and  thus  the  integration  of  Equations  (4.11)  and  (4.12)  becomes 

jyl  «■*•*  '-<*»  -  Wl 


where 


.at 

* 

i]  =  r  e’k°r**  dz  =  . 


d. 

ifr,  =  0 

sin(fcor,^) 

otherwise 

dy 

if  rw  =  0 

»a(fcory^) 

otherwise 

i*. 

3 

II 

*1  +  "j’ 

exp<jfcor,«i)  -  exp(-jtor,Tf-) 

jkorM 


otherwise. 


Due  to  the  product  en,_  P‘j_(0,d„j)  appearing  in  Equation  (4.11),  the  integrated 
electric  field  does  include  the  electric  field  resulting  from  the  open-current  end 
charge.  However,  since  the  product  enj_  P*j_(0,d„j)  occurs  outside  of  the  spatial 
integration  of  Equation  (4.11),  the  electric  field  of  the  open-current  end  charge  can 
be  easily  isolated  and/or  removed  using  the  technique  presented  in  Section  A.2.  This 
will  not  affect  the  values  of  /,,  Jy  or  i}. 
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4.2  =  0,  Region  (III) 

In  a  similar  manner,  1^^^  and  can  be  expressed 

gj*0*nj*+ 


tfm)  =  E  E  ^r1  ^+(o.<W  Xrn  *  (“*> 

fy >m  «**"♦*>.  (4.14) 


The  limits  of  integration  for  Region  (III)  are 


yni  _  j  _£  < 


_ii  <  .  <  ii 

a  -  *  -  a 

<  t 


a  s  y  s 

*3  <  *  <  £ 

and  thus  the  integration  of  Equations  (4.13)  and  (4.14)  becomes 

7 

/  e-i*0*  *+  _  /  /  rill 

yvme  aw  —  igif/ig 

where  7X  and  Iy  are  the  same  as  above  and 


-S3 


if  r,  =  0 
otherwise. 


Note  that  a  similar  statement  can  be  made  concerning  the  electric  field  resulting 
from  the  open-current  end  charge  for  the  integrated  Region  (III)  electric  field  as  was 
made  above,  for  the  Region  (I)  electric  field. 

4.3  ky,  =  0,  Region  (II) 

In  Region  (II),  the  spatial  variation  of  the  fields  consists  of  both  the  exponential 
term  e~'*°**f*,  and  also  the  pattern  factor  P*J±(o,6)  due  to  the  z  dependence  of  its 
integration  limits,  given  in  Section  A.l.  Thus,  and  £/n)  can  be  expressed  as 

T*(n)  _  ^  W 

^  24d,  trir 

—  J  fv n  *nj-Kj-(ani-Xj-)e-ikoR'*-<h  + 
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e*»*-i  f*  f  „  B  . 

+  — — - Jvu  «»j+  +  * 

< 


(4.15) 


T®(^)  _  1  y*  Y' 


-  r<  (*.,  x  »-)  /vII  ^J-(<^-Ai-)«-i*,I'-*>+ 

+  "  “  (^"J  x  f+)  ^#+(®iV+»Vf+)e  Afj  •  (^*W) 

The  electric  field  resulting  from  the  open-current  end  charge  is  embedded  in  the 
product  ^nj±  P^±(anj±%bnj±)‘  Thus  the  vector  e„j±  appears  inside  the  integral  of 
Equation  (4.15)  so  that  the  charge  electric  field  can  be  isolated  and/or  removed.  The 
limits  of  integration  for  Region  (II)  are 


<  »  <  ^ 

<  v  <  t 

<  i  <  *3, 


and  therefore  the  integrations  cl  Equations  (4.15)  and  (4.16)  can  be  written  as 

in  *"rt  *>  =  U. i3n> 

JvU  <fc  =  /.  V.±(II) 

where  /,  and  /„  are  the  same  as  above,  and 


ifi11’  =  P  *.jiPii±(.a«±,Ki±)t**‘r«dz 

/S(II> = r 

<1*1 
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For  the  evaluation  of  these  integrals,  the  nn&tion  of  Appendix  A  will  be  utilised, 
as  follows: 


Figure  A. 1(b)  Geometry 

Figure  A. 1(c)  Geometry 

-®n|+(Vj+) 

—  Bnj+i^j  +  finj  Z) 

•®nj+(^"J+) 

“  Q* nj+ 

■®ni+(a"i+) 

=  c;ni+ 

•®nj+(°nj+) 

—  Bnj+(<Znj  +  finj  *) 

Qnj+(bnj+) 

=  Qnj+i^j  +  finj  *) 

Qnj+(bnj+) 

=  Qnj+ 

Qnj+  (a»»i+  ) 

=  Q‘*h  =  o 

Qnj+{anj+) 

=  Qnj+{anj  +  finjz) 

^nj+(^j+) 

=  ^nj+(af»j  +  finj  Z) 

Anj+{bnj+) 

=  ^2n>+ 

^»i+(®"i+) 

=  Cinj+ 

i4*j+ (<*»»*+) 

=  ■^nj+(Q[t>j  +  @nj  «) 

B'i-ibnj-) 

=  C{nj. 

^i-(6ni-) 

=  B'^Onj  +  finj  Z) 

—  +  finj  z) 

Bn j-(°"j-) 

=  C[ni. 

Q'nj-fr nj-) 

=  Q'nj- 

Q'nj-Qnj-) 

—  Qnj—  (a«J  +  finjz) 

Qnj—  (anj-  ) 

=  Qnj—{°tnj  finj  z) 

Qnj-  (^nj—  ) 

=  Qnj-  =  0 

■^nj-(^nj-) 

=  c;nj_ 

•^nj-OW-) 

=  ■^nj-(ani  +  finj  z) 

^ni-(°»U-) 

—  Anj-(<*nj  +  finj  Z) 

■^»i-(°"i-) 

=  @2nj  -* 

(4-1 

The  terms  C'„i± 

and  Cjnj±  are  constants  with  respect  to 

z  which  result  from  t 

endpoints  of  segment  nj.  Similarly,  the  terms  Q„j±  are  the  charge  terms  associated 
with  the  open-current  at  the  endpoints  of  segment  nj,  and  are  independent  of  z.  The 
charge  terms  Qhj±(anj  +  finj  z)  must  be  removed  from  the  integral  for  the  answer 
to  be  right,  but  I’m  not  sure  why  this  is  so.  Using  the  definition  of  enj±  from 
Equation  (2.54),  and  Pnj±(&nj±,bnj±)  from  Appendix  A,  along  with  the  notation  of 
Equation  (4.17) 

ifi11’  =  +  «*]/£  + 

+  CJUSfjf* ->'*-<(** •**)]&+  <4-18> 

+  (Cf„,±  +  «/±)  +  i.,  C;„l±  ]  /3 

/.±(II)  =  cr±  S  [>•  (*.,  •  U)l}}±  -  llJ±  -  (4.19) 
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where 


where 


and 


5  = 


+1  if  Figure  A. 1(b)  geometry  and  +  integration 
or  Figure  A.  1(c)  geometry  and  —  integration, 
—1  if  Figure  A. 1(b)  geometry  and  —  integration 
or  Figure  A.l(c)  geometry  and  +  integration, 

/£.  =  r  e*jkor‘x  eiM-ni+Au'KA-j  **)  c08  [jto(crn>  +  0njz)\  dz 


_  gjko&nj(&nj't±) 


gj*0‘ 


-/+  + 


g-jtoanj 


I. 


l}}±  =  r  e*’kcr,x  sin  [fco(an;  +  A,;*)]  dz 

Jx  1 


gJ^OQfij  --jfcoanj 

/+  — — /- 


2j 


2 j 

**»  °nj*  ~F  yy  anjy 
Qniz 

if  =  -0, 
otherwise 


wnj  —  -F*1!*  *F  (dnj  *  f±)  — 

I±  =  [**  fa 

J*  i 


*3  “*1 


'»V 


/2  =  r  e*ikor”  dz  = 
Jx  l 


«3  —  Zi  if  r,  =  0 

safeihrsalzgefeattZial  otherwise. 


(4.20) 


(4.21) 


If  it  is  desired  to  remove  the  electric  held  resulting  from  the  open-current  end  charge, 
then  Q^j±  needs  to  be  removed  from  Equation  (4.18). 

4.4  kw  >  0  and  kw  <  0 

The  only  spatial  variation  of  the  fields  of  the  planar  arrays  where  k*  >  0  is  the 
exponential  term  e--**0®'*-  inside  the  double  spectral  summation  of  plane  waves. 
Thus  I®  (fc«,  >  0)  and  >  0)  can  be  expressed  as 

>0)  =  2i£££  (l -«-*-)  r,  X 
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X  *«j-  Kj-(0,<lnj)  Jv9  e*i*aE  #-  dv 
tH/i  1  e-*-  \  e****f- 


(4.M) 


x  (i^  x  *>•  (4.23) 

A  similar  statement  can  be  made  about  the  fields  of  the  planar  arrays  where  A*  <  0, 
and  thus  <  0)  and  <  0)  can  be  expressed  as 


^*-<0)=2&??(T^)^X 

I.''*'™**’ 

tH/i  «\  1  v-v-/  e+^+  \  ei*0,l^t+ 

x  (i^  x  f+)^+(0.4«)  /yt  «'AM*  <to. 


(4.24) 


(4.2S) 


The  integration  limits  are  over  the  entire  Q  =  0  lattice  cell  and  are  symmetric,  i.e., 

\  -*  <  *  <  £ 


V°  -  -**■<  y  <  % 

—4*  <  z  <  4m. 

3  -  Z  -  2  * 

The  integrations  of  Equations  (4.22)-{4.25)  can  be  written  as 

/,  «-*"**- 

where  J*  and  /v  are  the  same  as  above  and 

I,  =  [*  e*jk°r'’  dz=  d*  if  r,  =  0 

''“IT  sin(fcor,^)  otherwise. 

The  electric  field  of  open-current  end  charges  can  be  isolated  and/or  removed  by 
operating  on  the  enj±  P*J±(0,dnj)  product  as  outlined  in  Section  A.2. 
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Chapter  5 

Numerical  Results 

5.1  A  3D  Array  of  Short  PEC  Dipoles 

The  first  set  of  data  illustrate  the  convergence  of  the  effective  permittivity  and  the 
current  shape  of  a  simple  artificial  dielectric  with  respect  to  the  number  of  MM 
expansion  functions.  The  dhta  also  illustrate  the  ability  of  the  PMM  solution  to 
account  for  the  mutual  coupling  effects  between  wire  objects  in  the  artificial  dielectric. 
Finally,  it  is  shown  that  the  effective  wavenumber  can  vary  greatly  with  the  direction 
of  propagation,  whereas  the  effective  permittivity  is  typically  independent  of  the 
direction  of  propagation,  provided  the  current  shape  does  not  change  very  much. 

As  shown  in  the  insert  to  Figure  5.1,  the  geometry  of  the  artificial  dielectric 
consists  of  a  3D  periodic  array  of  short  perfect  electric  conducting  (PEC)  dipoles  in 
a  host  medium  of  free  space.  The  dipoles  have  radius  a  =  O.OOlAo  and  length  2 h  = 
0.2Ao,  and  are  oriented  parallel  to  the  z-axis.  They  are  arranged  in  a  lattice  structure 
where  dx  =  2h  +  d  and  dy  =  dx  =  d.  For  this  uniaxial  structure,  the  only  non-unity 
diagonal  element  of  the  permittivity  tensor  is  the  solution  will  have  an  x 

polarized  electric  current  and  electric  field.  The  plane  wave  is  propagating  in  the  fi*  = 
-(y+z)/v/2  direction.  The  computed  relative  effective  permittivity  is  plotted  for 
lattice  spacing  d  ranging  up  to  O.lAo  with  N  —  1,7  and  15  MM  expansion  functions 
for  the  current  on  each  wire  dipole.  The  N  =  1  curve  corresponds  identically  with  the 
results  obtained  by  Blanchard  [21,  Ch.  5].  Also,  an  N  =  1  static  approximation  for 
the  relative  effective  permittivity  is  given.  This  static  approximation  is  based  upon 
methods  presented  by  Collin  [20,  Ch.  12]. 
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Figure  5.1:  Relative  effective  permittivity  versus  lattice  spacing  for  a  3D  array  of 
perfectly  conducting  short  dipoles. 

Figure  5.2  shows  the  shape  of  the  normalized  current  induced  on  the  same  3D 
array  of  dipoles,  at  the  specific  value  of  d  =  0.02Ao,  for  N  =  1,7  and  15.  Note  that 
the  PMM  solution  predicts  a  rounded  off  pulse  shape  for  the  current,  and  several  MM 
PWS  expansion  functions  are  needed  to  accurately  model  the  current.  This  illustrates 
that  the  moment  method  solution  accounts  for  a  change  in  shape  in  the  dipole  current 
caused  by  mutual  coupling  effects  in  the  3D  array.  In  contrast,  typically,  static 
solutions  employ  a  polarizability  of  the  objects.  This  polarizability  is  a  characteristic 
parameter  of  the  objects,  and  is  not  influenced  by  the  3D  array  [20,  Sec.  12.1]. 
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I  =  Norm.  Current 


Figure  5.2:  The  normalized  current  induced  on  the  dipoles  in  a  3D  array  of  perfectly 
conducting  short  dipoles. 


57 


Now  consider  the  3D  array  of  PEC  dipoles  in  free  space  shown  in  the  insert  to 
Figure  5.3.  As  before,  the  dipoles  have  radius  a  =  0.001Ao  and  length  2k  =  0.2A<>. 
This  time  they  are  arranged  in  a  fixed  lattice  where  dx  =  0.25Ao  and  dy  =  dt  =  0.05Ao. 
Again,  the  artificial  dielectric  is  uniaxial,  and  e^x  is  the  only  non-unity  diagonal 
permittivity  tensor  component.  The  direction  of  propagation  U*  varies  with  the  angle 
8  from  0  to  80°,  measured  from  the  z-axis,  as  indicated  in  the  insert  to  Figure  5.3. 
The  data  in  Figures  5.3  and  5.4  were  computed  with  7  MM  expansion  functions. 

Figure  5.3  shows  the  relative  effective  permittivity  t"  (computed  by  both  the  po¬ 
larization  method  of  Equations  (2.32)  or  (2.38),  and  the  Maxwell’s  equations  method 
of  Equation  (2.37))  and  the  normalized  effective  wavenumber  (ke/ko)  as  a  function  of 
the  propagation  angle  8.  These  data  illustrate  that,  in  general,  the  effective  permit¬ 
tivity  tensor  components  will  be  almost  constant,  whereas  the  effective  wavenumber 
can  vary  quite  noticeably  with /different  propagation  directions.  This  is  exactly  true 
for  real  anisotropic  media,  and  it  holds  reasonably  well  for  many  typical  artificial 
anisotropic  media.  See  Appendix  G  for  a  discussion  of  the  wavenumber  predicted  by 
the  ellipsoid  of  wave  normals.  Figure  5.4  shows  the  current  shape  on  the  center  dipole 
for  various  angles  8.  Note  that  the  magnitude  of  the  current  shape  is  essentially  the 
same  for  different  propagation  direction  directions,  whereas  the  phase  of  the  current 
shape  exhibits  a  slight  odd  symmetry  about  the  center  of  the  dipoles. 

5.2  Dispersion  of  a  3D  Array  of  Dipoles 

The  data  in  Figure  5.5  illustrate  the  dispersion  characteristics  of  an  artificial  dielectric 
composed  of  PEC  dipoles.  As  shown  in  the  insert  to  Figure  5.5,  the  dipoles  are  of 
length  2 h  =  0.6cm,  radius  a  =  0.01cm,  and  are  arranged  in  a  3D  lattice  array 
with  spadngs  dx  =  0.7cm  and  dy  =  dx  =  0.1cm.  The  dipoles  are  parallel  to  the 
z-axis,  and  therefore  the  electric  current  and  electric  field  are  i  polarized.  This 
artificial  dielectric  is  uniaxial,  and  the  only  non-unity  diagonal  component  of  the 
permittivity  tensor  is  e^.  The  direction  of  propagation  is  along  the  z-axis.  The 
relative  effective  permittivity  f%x  is  plotted  versus  frequency  for  N  =  1,5  and  11  MM 
expansion  functions  for  the  current  on  each  dipole.  The  frequency  varies  up  to  24 
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Figure  5.3:  The  normalized  effective  wavenumber  (ke/ko)  and  the  relative  effective 
permittivity  versus  propagation  angle  9 ,  for  an  array  of  perfectly  conducting  short 
dipoles. 


Current  Phase  Current  Magnitude 


Figure  5.4:  Magnitude  and  phase  of  the  current  on  the  center  dipole  for  several 
different  propagation  angles  0,  for  an  array  of  perfectly  conducting  short  dipoles. 


60 


GHz,  corresponding  to  a  dipole  length  of  0.48A<j.  At  low  frequencies,  the  effective 
permittivity  approaches  a  constant  value,  and  the  N  =  1  solution  approaches  the 
static  approximation  of  Collin  [20].  As  the  frequency  increases,  the  scattered  field  of 
each  dipole  increases,  and  thus  the  effective  permittivity  increases,  especially  as  the 
half  wavelength  resonance  of  the  dipoles  is  approached. 

5.3  Array  of  Dipoles  Oriented  at  Angle 

The  data  in  this  section  show  an  example  of  a  simple  artificial  dielectric  that  is  not 
uniaxial,  i.e.,  the  effective  permittivity  tensor  has  non-zero  off-diagonal  components. 
It  is  shown  that  the  two  eigenfunction  solutions  (one  being  a  solution  with  a  free 
space  root  &*  =  fco)  are  needed  to  solve  for  the  permittivity  tensor.  As  shown  in  the 
insert  sketch  to  Figure  5.6,  the  artificial  dielectric  consists  of  straight  PEC  dipoles  of 
length  2 h  =  0.2Ao  and  radius  a  =  0.001  Ao, .oriented  in  the  direction  of  the  angle  <f>. 
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The  dipoles  are  arranged  in  a  3D  lattice  array  with  spacing!  dx  =  dy  =  0.2lAo  and 
dM  =  O.OIAq.  The  direction  of  propagation  is  along  the  z-axis. 

Figure  5.6  shows  the  normalized  effective  wavenumber  (he/ho)  as  a  function  of  the 
orientation  angle  for  N  =  7  and  15  MM  expansion  functions.  (Recall  that  another 
free  space  root  solution  exists,  too.)  Note  that  for  the  indicated  effective  wavenumber, 
the  electric  dipole  moment  P®  and  the  average  eigenfunction  electric  field  E®  are 
oriented  parallel  to  the  dipoles.  The  magnetic  dipole  moment  M®  vanishes  and  the 
average  eigenfunction  magnetic  field  H?  is  oriented  perpendicular  to  the  dipoles  in 
the  xy-plane.  Corresponding  to  the  free  space  root  he  =  A^>,  the  dipole  moments 
(P®,  M®)  vanish  and  the  eigenfunction  fields  may  be  written  as 

Ej  =  xrfo  sin  <f>  —  ytjo  cos  4> 

jH®  =  x  cos  +  y  sin  0  (5.1) 

where  tjq  is  the  characteristic  impedance  of  free  space.  Using  the  current  moment  and 
eigenfunction  fields  above  in  Equations  (2.32)  or  (2.37),  the  polarization  method  and 
the  Maxwell’s  equations  method  yield  virtually  identical  results  for  the  permittivity 
tensor  components,  shown  in  Figure  5.7.  Note  that  the  artificial  dielectric  degenerates 
to  a  uniaxial  dielectric  when  the  dipoles  are  oriented  along  a  principle  axis,  i.e.,  when 
4>  =  0  or  <f>  =  90°. 

5.4  Resistive  Loaded  Dipoles 

This  data  in  this  section  show  the  effect  of  resistive  loading  on  a  3D  array  of  dipoles. 
As  shown  in  the  insert  to  Figure  5.8,  the  geometry  consists  of  PEC  dipoles  of  length 
2 h  =  0.2Ao  and  radius  a  =  O.OOlAo,  arranged  in  a  3D  lattice  with  spadngs  dx  = 
0.23Ao  and  dy  =  dx  =  0.03Ao.  A  purely  resistive  load  Rl  is  located  at  the  center 
of  each  dipole.  Propagation  is  in  the  iifc  =  z  direction,  and  polarization  is  in  the 
x  direction.  Figure  5.8  shows  the  relative  effective  permittivity  and  effective 
loss  tangent  tan  Sesx  of  the  artificial  medium  versus  load  resistance  for  10Q  <  Rl  < 
100KH,  for  solutions  with  N  —  7  and  N  =  13  MM  expansion  functions.  Note  that 
at  low  resistive  loading  (Rl  <  200ft)  the  effective  permittivity  is  close  to  that  for 
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Figure  5.6:  The  normalized  effective  wavenumber  (he/ho)  versus  orientation  angle 
for  an  array  of  perfectly  conducting  short  dipoles. 
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Relative  Permittivity  Tensor  Components 


Figure  5.7:  The  effective  permittivity  tensor  components  versus  orientation  angle  <f> 
for  an  array  of  perfectly  conducting  short  dipoles. 
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rl  =  Load  Resistance  in  Ohms  ( Q ) 


Figure  5.8:  Relative  effective  permittivity  and  loss  tangent  for  a  3D  array  of  resistive 
loaded  PEC  dipoles. 

PEC  dipoles  of  length  2h,  and  at  high  resistive  loading  (Rl  >  10K ft)  the  effective 
permittivity  is  close  to  that  for  two  disconnected  PEC  dipoles,  each  of  length  h.  The 
maximum  loss  of  the  artificial  medium  occurs  near  Rl  =  1.5KQ,  corresponding  to 
the  maximum  PRl  loss  in  the  load  resistance  of  the  dipoles. 

Figure  5.9  shows  the  magnitude  and  phase  of  the  current  shape  on  the  center  dipole 
for  N  =  13  MM  expansion  functions  and  for  load  resistances  of  Rl  =  100ft,  lKfl  and 
lOKft.  Note  that  as  the  load  resistance  increases,  the  current  at  the  terminals  of  the 
dipole  (at  Rl)  decreases,  until  for  large  Rl,  the  dipole  is  essentially  disconnected  at 
its  center,  effectively  forming  two  dipoles. 
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Figure  5.9:  Magnitude  and  phase  of  the  current  on  the  center  dipole  at  Ri  =  100(1, 
lKfl  and  10KQ,  for  a  3D  array  of  resistive  loaded  PEC  dipoles. 
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5.5  Lossy  Dielectric  Dipoles 

The  data  in  this  section  show  the  effects  of  dielectric  loss  in  a  3D  array  of  dipoles.  For 
simplicity  of  dimensions,  the  data  are  computed  at  300  MHz  with  a  free  space  host 
medium  so  that  Ao  =1  meter.  As  shown  in  the  insert  to  Figure  5.10,  the  geometry 
consists  of  lossy  dielectric  dipoles  of  length  2 h  =  0.2Ao  and  radius  a  =  O.OOlAo, 
arranged  in  a  3D  lattice  with  spacings  dx  =  0.23Ao  and  dy  =  dM  =  0.03Ao.  The 
dipoles  have  relative  dielectric  constant  cJr  =  1  and  loss  tangent  tan£i.  Propagation 
is  in  the  u*  =  z  direction,  and  polarization  is  in  the  x  direction.  Figure  5.10  shows 
the  relative  effective  permittivity  e^x  and  effective  loss  tangent  tan  8XX  of  the  artificial 
medium  versus  dipole  loss  tangent  for  10  <  tan  8\  <  100, 000  and  for  solutions  with 
N  =  9  and  N  =  15  MM  expansion  functions.  Note  that  a  monopole  expansion 
function,  with  its  associated  open-current  end  charge  contribution,  was  included  at 

J 

each  end  of  the  dielectric  dipoles. 

Figure  5  .11  shows  the  magnitude  and  phase  of  the  current  shape  on  the  center 
dipole  at  N  =  15  MM  expansion  functions  for  dipole  loss  tangent  values  of  tan  Si  = 
50,  500  and  5000.  Note  that  for  a  low  dipole  loss  tangent,  the  current  is  fairly  uniform 
across  the  dipole.  As  the  dipole  loss  tangent  increases,  the  current  shape  approaches 
that  of  a  PEC  dipole. 

5.6  Square  PEC  Loops 

The  data  in  this  section  involves  an  artificial  medium  that  is  slightly  magnetic.  It 
is  shown  that  the  eigenfunction  current  shape  on  the  wire  objects  can  be  a  strong 
function  of  the  direction  of  propagation.  As  shown  in  the  insert  to  Figure  5.12,  the 
geometry  of  the  artificial  medium  consists  of  a  3D  periodic  array  of  PEC  square  loops 
arranged  in  a  host  medium  of  free  space.  The  loop  wires  have  radius  a  =  O.OOlAo  and 
the  sides  are  of  length  l  =  O.lAo.  The  loops  are  arranged  perpendicular  to  the  z-axis 
in  a  lattice  structure  where  dx  =  dy  =  0.12Ao  and  dx  =  0.02Ao.  Propagation  is  in  the 
U*  direction,  always  in  the  the  yz-plane,  measured  by  the  angle  6  from  the  z-axis.  6 
varies  from  0°  (broadside  to  the  the  loops)  to  90°  (edge-on  to  the  loops).  For  all  the 
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Figure  5.11:  Magnitude  and  phase  of  the  current  on  the  center  dipole  at  tan  Si  =  50, 
500  and  5000,  for  a  3D  array  of  lossy  dielectric  dipoles. 
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data  shown  in  this  section,  there  are  8  MM  expansion  functions  distributed  equally 
around  the  square  loop. 

Figure  5.12  shows  the  normalized  effective  wavenumber  (fce/ho)  corresponding  to 
each  of  the  two  roots  (see  Section  2.3.1)  as  a  function  of  the  propagation  angle  0.  The 
roots  were  computed  by  two  methods.  The  first  method  determines  the  roots  by  the 
PMM  solution.  The  second  method  employs  a  technique  from  crystal  optics  known  as 
the  ellipsoid  of  wave  normals  (see  Appendix  G)  [23,  Ch.  4], [24,  Ch.  14).  The  ellipsoid 
of  wave  normals  uses  the  roots  computed  by  the  PMM  for  0  =  0°  and  0  =  90°. 
Figure  5.13  shows  the  magnitude  of  the  determinant  of  the  impedance  matrix  \Z\ 
versus  the  normalized  effective  wavenumber  for  propagation  angles  0  =  0, 15, 30, 60 
and  90°.  The  roots  are  indicated  by  the  sharp  minima  of  \Z\.  Note  the  occurrence 
of  a  double  root  due  to  symmetry  considerations  only  at  0  =  0°.  Also,  note  that 
the  0  =  90°  curve  only  indicates  one  root;  the  other  root  being  a  free  space  root  at 
fee  =  Jfeo. 

Figures  5.14  and  5.15  show  the  current  shape  on  the  center  loop  for  both  roots 
1  and  2,  respectively,  for  several  propagation  angles  9.  The  sketch  at  the  top  of 
Figure  5.15  shows  how  the  current  is  plotted  in  relation  to  the  geometry  of  the  loop. 
The  phase  of  the  current  shape  corresponding  to  root  2  is  not  plotted  because  it  is 
0°  for  positive  z  and  180°  for  negative  z.  Note  that  the  current  shape  corresponding 
to  root  1  changes  with  propagation  angle  0,  whereas  the  current  shape  corresponding 
to  root  2  does  not  change  with  propagation  angle.  Also,  note  that  at  propagation 
angle  0  =  0°  the  two  current  shapes  are  orthogonal  to  each  other,  corresponding  to 
the  double  root  from  symmetry.  The  current  shape  for  root  1  has  net  electric  dipole 
moment  P°  oriented  in  the  z-direction,  whereas  the  current  shape  for  root  2  has  P° 
oriented  in  the  y-direction. 

Figures  5.16  and  5.17  show  the  effective  relative  permittivity  and  permeability 
tensor  components,  respectively,  computed  as  a  function  of  the  propagation  angle  0. 
Observe  that  this  artificial  medium  is  uniaxial.  The  permittivity  tensor  components 
were  computed  using  both  the  polarization  method  of  Equations  (2.32)  or  (2.38), 
and  the  Maxwell’s  equations  method  of  Equation  (2.37).  There  is  good  agreement 
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Figure  5.12:  Normalized  effective  wavenumber  (fce/fco)  for  both  roots  versus  propaga¬ 
tion  angle  0,  for  an  array  of  square  PEC  loops. 
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Figure  5.13:  Magnitude  of  \Z\  versus  normalized  effective  wavenumber  at  6 
0, 15,30,60  and  90°,  for  an  array  of  square  PEC  loops. 


1.2 

a  1.0 

i  0.8 

o>  0.6 

s 

0.4 

Jn 

9  0.2 

u 

0.0 

0 

d>  -20 

a> 

8  -40 

S  "60 

5  -80 

04 

.  -100 
M 

£  -120 

°  -140 


- —  =  Distance  Around  Loop 
A0 


Figure  5.14:  The  root  1  current  mode  shape  at  9  =  0, 10, 20, 30, . . .  ,90°,  for  an  array 
of  square  PEC  loops. 
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Figure  5.15:  The  root  2  current  mode  shape  at  all  9 ,  for  an  array  of  square  PEC 


between  the  two  methods  for  and  e£,  but  not  for  Note  that  and 
are  solely  determined  from  the  root  2  solution,  whose  current  shape  does  not  change 
with  propagation  angle  9.  On  the  other  hand,  is  determined  solely  from  the  root 
1  solution,  where  the  current  shape  does  change  with  propagation  angle  9.  This 
discrepancy  between  the  two  solutions  for  e^.  has  not  been  resolved.  However,  it  is 
interesting  to  note  one  possible  major  difference  between  real  anisotropic  media  and 
artificial  media.  In  a  real  anisotropic  medium,  the  permittivity  tensor  components 
will  not  change  with  propagation  direction.  In  an  artificial  medium,  the  permittivity 
tensor  components  may  change  with  propagation  direction  since  the  eigenfunction 
currents  may  change. 

5.7  PEC  Wire  Crosses 

The  data  in  this  section  illustrate  the  effect  of  the  horizontal  cross  member  on  the 
dispersion  characteristics  for  an  array  of  PEC  wire  crosses.  As  shown  in  the  insert 
to  Figure  5.18,  the  wire  crosses  have  vertical  extent  from  y  =  —L/2  to  y  ==  +1/ 2, 
and  horizontal  extent  from  *  =  —L/A  to  as  =  +L/A,  with  L  =  5cm.  The  horizontal 
cross  member  is  located  at  y  =  +L/ 4,  and  the  wire  radius  is  a  =  1mm.  The  wire 
crosses  are  arranged  in  a  3D  lattice  with  dx  —  3.75cm,  dy  =  7.5cm  and  dt  =  1cm. 
This  artificial  dielectric  is  uniaxial  with  non-unity  values  for  both  and  e^.  The 
direction  of  propagation  is  along  the  z-axis.  There  are  7  vertical  expansion  functions, 
3  horizontal  expansion  functions,  and  1  horizontal  to  vertical  expansion  function,  for  a 
total  of  N  =  11.  Data  are  also  included  for  the  vertical  and  horizontal  members  alone 
(arrays  of  dipoles),  as  shown  in  the  inserts  for  the  lower  two  curves  of  Figure  5.19. 

Figure  5.18  shows  the  dispersion  variation  of  e^x  and  for  frequencies  up  to  3 
GHz,  corresponding  to  L  =  0.5Ao,  for  the  3D  array  of  wire  crosses.  Also,  shown  are 
the  dispersion  variations  of  for  the  array  of  horizontal  dipoles,  and  for  the  array 
of  vertical  dipoles.  Corresponding  to  the  horizontally  polarized  (x-polarized)  wave, 
note  that  e^.  is  identical  for  the  array  of  wire  crosses  and  the  array  of  horizontal 
dipoles.  The  horizontal  eigenfunction  current  is  the  same  for  both  cases  because 
the  vertical  member  is  symmetrically  located  with  respect  to  the  horizontal  member. 
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Relative  Permittivity  Tensor  Components 


Figure  5.16:  The  permittivity  tensor  components  versus  propagation  angle  8,  for  an 
array  of  square  PEC  loops. 
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Relative  Permeability  Tensor  Components 
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Figure  5.17:  The  permeability  tensor  components  versus  propagation  angle  0 ,  for  an 
array  of  square  PEC  loops. 
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Frequency  (GHz) 


Figure  5.18:  Dispersion  curve  for  a  3D  array  of  PEC  wire  crosses. 


However,  for  vertical  polarization,  ej  is  greater  for  the  array  of  wire  crosses  than  for 
the  array  of  vertical  dipoles.  The  horizontal  member  is  not  symmetrically  located 
with  respect  to  the  vertical  member,  allowing  for  a  greater  net  y -directed  current 
in  the  case  of  the  wire  crosses.  This  effect  can  also  be  observed  in  Figure  5.19, 
which  shows  the  magnitude  of  the  determinant  of  the  impedance  matrix  \Z\  versus 
the  normalized  effective  wavenumber  for  the  array  of  PEC  wire  crosses,  as  well  as 
for  the  arrays  of  vertical  and  horizontal  dipoles,  at  the  frequency  of  2  GHz.  Note 
that  the  horizontal  polarization  root  does  not  change  from  the  array  of  wire  crosses 
to  the  array  of  horizontal  dipoles,  whereas  the  vertical  polarization  root  occurs  at  a 
greater  effective  wavenumber  for  the  array  of  wire  crosses  than  for  the  array  of  vertical 
dipoles. 


5.8  Graphite-Epoxy  2D  Composite  Medium 


This  section  presents  the  analysis  of  a  modern  composite  material  consisting  of  very 
thin  graphite  fibers  embedded  in  an  epoxy  host  binding  material.  Figure  5.20  shows 
the  geometry  of  this  composite  material.  The  graphite  fibers  are  modeled  as  material 
wires  of  infinite  length  in  the  x  direction  with  radius  a  =  3.2pm,  spaced  in  a  square 
2D  lattice  with  dy  —  dt  =  7.5pm.  The  conductivity  of  the  graphite  fibers  is  71.4 
K(fl~1)/meter  and  the  permittivity  of  the  host  epoxy  material  is  cor  =  2.5.  As  shown 
in  Figure  5.20,  the  PMM  solution  uses  3  expansion  functions;  1  dipole  expansion 
function  and  two  monopole  expansion  functions  without  open-current  end  charges. 
(The  charges  cancel  between  adjacent  x  direction  cells.)  To  model  this  2D  material, 
the  x  direction  lattice  spacing  was  chosen  as  ds  =  75pm. 

Figure  5.21  shows  the  computed  dispersion  characteristics  of  the  composite  ma¬ 
terial.  Note  that  at  low  frequencies,  the  effective  conductivity  is  very  close  to  what 
results  from  using  a  simple  fill  factor  formula  based  on  a  ratio  of  the  area  occupied  by 
the  graphite  fibers  to  the  area  occupied  by  a  2D  lattice  cell,  i.e.,  for  low  frequencies, 
the  effective  conductivity  is  approximately  given  by 


rar 


(5.2) 
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Figure  5.19:  Magnitude  of  \Z\  versus  normalized  effective  wavenumber  for  an  array 
of  PEC  wire  crosses,  and  for  the  isolated  vertical  and  horizontal  dipole  members. 
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Figure  5.20:  The  geometry  of  the  graphite-epoxy  2D  composite  medium. 
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Applying  Equation  (5.2)  to  the  graphite-epoxy  composite  medium  results  in  ffj,  = 
40.8  K(fl_1)/meter,  agreeing  very  closely  with  the  results  of  Figure  5.21. 


5.9  Dielectric  Weave 

The  data  in  this  section  show  the  dispersion  characteristics  of  the  effective  permittivity 
of  a  dielectric  weave ,  a  geometry  which  has  current  that  flows  between  adjacent  lattice 
cells.  As  shown  in  the  insert  to  Figure  5.22,  the  geometry  consists  of  stacked  or 
layered  square  grids  of  dielectric  rods.  The  dielectric  rods  have  relative  permittivity 
cjr  =  10,  loss  tangent  tan  6\  —  0  or  1,  and  radius  a  =  2mm.  The  grid  dimensions  are 
ds  —  dy  —  L  —  5cm,  and  are  spaced  a  distance  of  dM  =  6mm  apart.  Propagation  is 
along  the  z-axis,  and  due  to  symmetry  considerations,  the  medium  is  uniaxial  with 
and  =  1.  All  data  was  computed  with  N  =  11  expansion  functions, 
including  4  monopole  expansioh  functions  to  enforce  continuity  of  current  between 
adjacent  lattice  cells.  Figure  5.22  shows  the  relative  effective  permittivity  e^.  = 
and  effective  loss  tangent  tan  6*,  =  tan  6^  of  the  artificial  medium  for  frequency 
varying  up  to  3  GHz  (corresponding  to  a  grid  size  of  L  =  0.5Ao)  for  dielectric  rod 
loss  tangent  values  of  tan  Si  =  0  and  1.  Note  that  the  relative  effective  permittivity 
and  effective  loss  tangent  are  almost  constant  across  the  given  frequency  range.  This 
is  due  to  the  fact  that  the  current  is  essentially  constant  in  the  dielectric  rods  from 
continuity  of  current  between  adjacent  lattice  cells. 
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Figure  5.21:  Dispersion  curves  for  a  composite  graphite-epoxy  material. 
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Chapter  6 

Computer  Program  ADWIRS 


This  chapter  presents  the  usage  of  the  computer  program  ADWIRS.  The  program 
ADWIRS  can  analyze  an  artificial  medium  composed  of  a  3D  periodic  array  of 
identical  arbitrarily-shaped  thin  conductive  or  dielectric  wire  objects  arranged  in  a 
homogeneous  host  medium.  >  ADWIRS  computes  such  parameters  as  the  effective 
wavenumber  corresponding  to  the  root,  the  shape  of  the  eigenfunction  currents 
in  the  wire  objects,  the  average  electric  and  magnetic  dipole  moments  per  unit  cell 
(P°,M°),  and  the  average  electric  and  magnetic  fields  per  unit  cell  (E°,H°). 

Concerning  the  effective  constitutive  parameters  of  the  artificial  medium,  AD¬ 
WIRS  only  computes  the  effective  permittivity  Ce  as  given  in  Equation  (2.41).  Note 
that  this  equation  for  evaluating  Ce  is  only  valid  when  polarization  and  direction  of 
propagation  are  along  one  of  the  principle  axes.  The  effective  permittivity  e«  can  be 
written  as 

=  =  (6.1) 

u> 

where  the  u  =  xx,yy  or  zz,  as  determined  by  the  principle  axis  of  polarization. 
ADWIRS  computes  the  three  values  e£,  tan  £?•  and  0^.  Often  Equation  (2.41)  does 
not  apply,  and/or  the  polarization  method  (Section  2.3.2)  and  Maxwell’s  equations 
method  (Section  2.3.3)  are  used  to  compute  the  effective  constitutive  parameters.  In 
this  case,  (P°,M°)  and  (E°,H°)  must  be  used  in  another  program  to  be  written  by 
the  user,  which  solves  for  ee  and  jie  using  the  polarization  method  or  the  Maxwell’s 
equations  method. 
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C  READ  1: 

Read  (8,*)  ConTol,Derr, Ipv, Iwrz, Iwrc, Iqz, Isum,Nmaxz,ltaaxe 
C  READ  2: 

Read  (8,*)  Iswp,Nswpts , Parami , Dparam 
C  READ  3: 

Read  ( 8 , * )  PGHz , Eps  Or , TandO , Epserl , Tandel , Epser2 , Tande2 
C  READ  4: 

Read  (8, *)  Np,Ns,Nload 
C  READ  5: 

Do  N  =  l,Np 

Read  (8,*)  X<N) ,Y(N) ,Z(N) 

End  Do 
C  READ  6: 

Do  N  =  l,Ns 

Read  (8,*)  Ia{N) , Ib(N) , Iloss(N) , Radwire(N) , Epsrwire(N) .Loss 
<  FORTRAN  statements  to  process  Loss  > 

End  Do 
C  READ  7: 

Do  N  =  l.Nload 

Read  (8, *)  Iload(N) , Zload(N) 

End  Do 
C  READ  8: 

V(3)  =  0.0 

Read  (8,  *)  V(1),V(2) 

C  READ  9: 

Read  (8,*)  W(l) ,W(2) ,W(3) 

C  READ  10: 

Read  (8, *)  D(l) ,D(2) ,D(3) 

C  READ  11: 

Read  (8,*)  Uk(l) ,Uk(2) ,Uk(3) 


Figure  6.1:  The  ADWIRS  program  READ  statements. 


6.1  Inputs  To  ADWIRS 


The  inputs  to  the  computer  program  ADWIRS  are  read  via  11  unformatted  FOR¬ 
TRAN  77  READ  statements  from  an  input  file  on  logical  unit  8.  These  READ 
statements  are  shown  in  Figure  6.1.  The  values  defined  in  the  READ  statements  are 
explained  in  the  sections  that  follow,  along  with  the  type  of  FORTRAN  variable  the 
input  value  must  be  in  the  input  file  (Integer,  Real  or  Complex.) 


6.1.1  READ  1:  Run  Control  Parameters 


ConTol  =  Convergence  test  tolerance  for  the  spectral  evaluation  of  the  impedance 
matrix  term  contributions  Z*,  given  in  Section  3.3,  and  for  the  average  eigen¬ 
function  fields  (E°,  H®),  given  in  Section  4.  Note  that  all  spectral  self  and  mu¬ 
tual  impedance  terms  Z£n  are  evaluated  to  within  a  tolerance  of  ConTol.  Also, 
the  average  eigenfunction  fields  (E°,H®)  have  converged  to  within  a  tolerance 
of  ConTol.  If  ConTol  is  set  greater  than  0.1,  then  ConTol  is  automatically 
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set  equal  to  0.1,  so  there  is  never  worse  than  10%  accuracy.  Typically  ConTol 
=  0.01,  corresponding  to  1%  accuracy.  (Real) 

Derr  =  Convergence  test  tolerance  for  the  evaluation  of  the  root  he,  i.e.,  the  root 
search  for  he  continues  until  he  is  determined  to  within  a  tolerance  of  Derr.  If 
Derr  is  set  greater  than  0.1,  then  Derr  is  automatically  set  equal  to  0.1,  so 
there  is  never  worse  than  10%  accuracy.  Typically  Derr  =  0.01,  corresponding 
to  1%  accuracy.  (Real) 

Ipv  =  Indicator  for  a  2D  planar  lattice  or  a  3D  volume  lattice  impedance  calculation. 
A  2D  planar  lattice  includes  only  the  ky,  =  0  planar  array  of  wire  objects.  A 
3D  volume  lattice  includes  all  the  planar  arrays  indexed  by  K,  (h*,  =  0,h«,  >  0 
and  kw  <  0).  (See  Section  3.3.)  (Integer) 

=  1  implies  only  compete  the  2D  planar  lattice  impedances  from  the  h*,  =  0 
plane  of  wire  objects. 

=  2  implies  compute  the  3D  volume  lattice  impedances  from  all  planar  array 
indexed  by  ky,  (kw  =  Ojh^  >  0  and  ky,  <  0). 

Note  that  for  propagation  through  a  3D  array  of  wire  objects,  set  Ipv  =  2. 

Iwrz  =  Indicator  for  printing  the  impedance  matrix.  (Integer) 

=  0  implies  do  NOT  print  the  impedance  matrix. 

=  1  implies  print  the  impedance  matrix. 

Note  that  setting  Iwrz  =  1  will  result  in  all  the  different  contributions  to  the 
impedance  matrix  being  printed  at  each  value  of  he  during  the  root  search 
iteration  procedure. 

Iwrc  =  Indicator  for  printing  the  current  solution.  (Integer) 

=  0  implies  do  NOT  print  the  current  solution. 

=  1  implies  print  the  current  solution. 

Note  that  setting  Iwrc  =  1  will  result  in  the  current  solution  being  printed  for 
all  possible  choices  of  current  coefficients  /„  set  to  unity  (see  Section  2.2.2). 
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Iqz  =  Indicator  for  including  open-current  end  charge  contributions  in  the  evaluation 
of  the  electric  fields  and  impedance  matrix  terms.  (Integer) 

=  0  implies  do  NOT  include  open-current  end  charge  contributions. 

=  1  implies  include  open-current  end  charge  contributions. 

Note  that  the  user  can  always  set  Ipv  =  1  to  include  open-current  end  charge 
contributions  and  get  theoretically  correct  results.  However,  including  open- 
current  end  charge  contributions  usually  makes  the  impedance  matrix  converge 
much  more  slowly,  and  thus  it  is  recommended  to  set  Ipv  =  0  when  possible.  In 
general,  the  user  can  set  Ipv  =  0  except  when:  1)  open-current  end  monopole 
basis  functions  exist  in  the  material  wire  scattering  object  (for  example  with 
dielectric  wires),  and  2)  these  open-current  end  monopole  basis  functions  do 
NOT  connect  with  a4joining  monopole  basis  function  from  an  adjacent  lattice 
cell  (for  example  with  a  dielectric  weave). 

Isum  =  Indicator  for  performing  an  approximate  spectral  summation  on  root  search 
values  beyond  the  initial  guess.  (Integer) 

=  0  implies  perform  exact  spectral  summation  on  all  root  search  values. 

—  1  implies  perform  approximate  spectral  summation  beyond  the  initial  root 
search  guess. 

The  user  can  always  set  Isum  =  0  to  obtain  a  solution.  However,  setting  Isum 
=  1  can  save  considerable  CPU  time,  without  much  loss  in  accuracy,  when  the 
initial  root  guess  is  dose  to  the  final  root.  Typically,  the  user  can  set  Ipv  =  1 
and  obtain  a  solution.  Then  the  user  should  check  that  the  final  root  is  within 
a  few  percent  the  initial  root  guess.  If  not,  the  user  should  use  the  final  root  of 
this  solution  as  the  initial  root  guess  of  a  new  solution,  and  then  compute  the 
new  solution. 

Nmaxz  =  Maximum  double  summation  spectral  index  limit  for  the  evaluation  of  the 
spectral  impedance  matrix  contributions  Z*,  i.e.,  n*  and  n*  in  Sections  3.3.1 
and  3.3.2  shall  not  exceed  Nmaxz,  even  in  the  event  that  all  the  spectral 
impedance  matrix  terms  Z^n  have  not  converged  to  ConTol  by  this  point. 
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Typically,  Nmaxz  =  50,  if  no  open-current  end  monopole  basis  functions  exist 
in  the  scattering  object  geometry.  If  open-current  end  monopole  basis  functions 
do  exist  in  the  scattering  object  geometry,  then  fields  converge  much  more 
slowly,  and  hence  Nmaxz  =  100  or  greater.  (Integer) 

Nmaxe  =  Maximum  double  summation  spectral  index  limit  for  the  evaluation  of 
the  average  eigenfunction  fields  (E°,H°),  i.e.,  n„  and  n„  in  Chapter  4  shall  not 
exceed  Nmaxz,  even  in  the  event  that  all  the  average  eigenfunction  fields  have 
not  converged  to  ConTol  by  this  point.  The  average  eigenfunction  fields  are 
much  easier  to  evaluate,  and  converge  more  quickly  than  the  spectral  impedance 
matrix  contributions,  and  hence  typically  Nmaxe  =  20.  (Integer) 

6.1.2  READ  2:  Parameter  Sweep  Inputs 

Iswp  =  Indicator  for  what  type  of  a  parameter  sweep  to  perform.  (Integer) 

=  —  1  implies  sweep  Jfee,  and  evaluate  impedance  matrix  determinant  only. 

=  0  implies  no  parameter  sweep,  i.e.,  evaluate  a  single  solution. 

=  1  implies  sweep  frequency  in  GHz. 

=  2  implies  sweep  du  (in  meters)  where  <£,,  and  d*,  remain  fixed. 

=  3  implies  sweep  d„  (in  meters)  where  du  and  d^  remain  fixed. 

=  4  implies  sweep  d«,  (in  meters)  where  du  and  dp  remain  fixed. 

=  5  implies  sweep  d„  =  d^  (in  meters)  where  d*,  remains  fixed. 

=  6  implies  sweep  du  =  du,  (in  meters)  where  d*  remains  fixed. 

=  7  implies  sweep  d„  =  d u,  (in  meters)  where  du  remains  fixed. 

=  8  implies  sweep  eor  =  relative  permittivity  of  the  host  medium. 

=  9  implies  sweep  tan  So  =  loss  tangent  of  the  host  medium. 

=  10  implies  sweep  a  =  radius  (in  meters)  of  every  wire  segment.  (All  wire 
segments  are  set  the  same.) 

=  11  implies  sweep  eir  =  relative  permittivity  of  every  wire  segment.  (All  wire 
segments  are  set  the  same.) 

=  12  implies  sweep  tan  6\  —  loss  tangent  of  every  wire  segment.  (All  wire 
segments  are  set  the  same.) 
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=  13  implies  sweep  <rx  =  conductivity  (in  Q-1/  meter)  of  every  wire  segment. 
(AH  wire  segments  are  set  the  same.) 

Nswpts  =  Number  of  points  in  the  parameter  sweep.  (Integer) 

Parami  =  Initial  value  of  the  parameter  being  varied.  (Real) 

Dparam  =  Increment  step  size  of  the  parameter  being  varied.  (Real) 

Note  that  if  an  impedance  matrix  determinant  sweep  is  performed  (Iswp  =  —  1) 

then  the  variables  Parami  and  Dparam  refer  to  normalised  wavenumber  values,  i.e., 

Jfee/Jfc&.  Also,  the  effective  wavenumber  K  must  be  a  purely  real  value. 

6.1.3  READ  3:  Frequency,  Host  Media  and  Initial  Root 
Guess 

/ 

FGHz  =  Frequency  in  GHz.  (Real) 

EpsOr  =  «br  =  relative  permittivity  of  the  host  medium.  (Real) 

TandO  =  tan  S0  =  loss  tangent  of  the  host  medium.  (Real) 

Epserl  =  Initial  guess  at  £„  =  relative  effective  permittivity,  corresponding  to  the 
initial  value  in  the  parameter  sweep  (if  Iswp  >  0)  or  at  the  single  computed 
value  (if  Iswp  =  0.)  (Real) 

Tfendel  =  Initial  guess  at  tan  Se  =  effective  loss  tangent,  corresponding  to  the  initial 
value  in  the  parameter  sweep  (if  Iswp  >  0)  or  at  the  single  computed  value  (if 
Iswp  =  0.)  (Real) 

Epser2  =  Initial  guess  at  e„.,  corresponding  to  the  second  value  in  the  parameter 
sweep  (if  Iswp  >  0.)  (Real) 

TUnde2  =  Initial  guess  at  tan  S„,  corresponding  to  the  second  value  in  the  parameter 
sweep  (if  Iswp  >  0.)  (Real) 
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Note  that  the  exact  complex  permittivity  of  the  host  medium  it 

eo  =  cor  e(l  -jtan^o) 

where  e  is  the  free  space  permittivity.  Also,  in  the  program  ADWIRS,  the  root  he 
is  related  to  the  effective  permittivity  as 

he  =  uty/fioCt  where  e*  =  e*,.  e(l  —  j  tan  Se). 

If  a  parameter  sweep  is  computed  (Iswp  >  0  in  READ  2)  then  the  initial  guess  at 
the  first  two  roots  are  defined  through  the  input  variables  Epserl,  Tandel,  Epser2 
and  TWnde2.  The  initial  guess  at  subsequent  roots  is  generated  automatically  within 
the  program  ADWIRS.  If  only  a  single  value  is  computed  (Iswp  =  0  in  READ  2) 
then  the  initial  guess  at  the  root  is  defined  by  Epserl  and  Tandel.  Finally,  if  an 
impedance  matrix  determinant  sweep  is  computed  (Iswp  =  —  1  in  READ  2)  then  no 
initial  root  guess  is  required,  but  Epserl,  Tandel,  Epser2  and  Tande2  must  be 
input  anyway  for  consistency. 

6.1.4  READ  4:  Number  of  Points,  Segments  and  Lumped 
Loads 

READ  statements  4,5  and  6  define  the  geometry  of  the  wire  objects.  The  method 
of  wire  geometry  input  used  by  ADWIRS  is  too  complicated  to  explain  here,  but 
has  been  used  before  and  is  well  documented.  The  reader  is  referred  to  [34],  [35] 
and/or  [36]  for  a  detailed  description  of  the  wire  geometry  input  method.  The  input 
parameters  are  briefly  explained  below. 

Np  =  Total  number  of  wire  paints  defining  the  scattering  object.  (Integer) 

Ns  =  Total  number  of  wire  segments  defining  the  scattering  object.  (Integer) 
Nload  as  Total  number  of  lumped  loads  defining  the  wire  object.  (Integer) 

6.1.5  READ  5:  Wire  Point  Coordinates 

X(N)  =  e-coordinate  (in  meters)  of  wire  point  N,  for  N  =  1,2,. ..,Np.  (Real) 
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Y(N)  =  y-coordinate  (in  meten)  of  wire  point  N,  for  N  =  l,2,...,Np.  (Real) 

Z(N)  =  z -coordinate  (in  meten)  of  wire  point  N,  for  N  =  l,2,...,Np.  (Real) 

6.1.6  READ  6:  Wire  Segments 

Ia(N)  =  Endpoint  A  of  wire  segment  N,  for  N  =  l,2,...,Ns.  (Integer) 

Ib(N)  =  Endpoint  B  of  wire  segment  N,  for  N  =  l,2,...,Ns.  (Integer) 

Hoaa(N)  =  Indicator  for  the  material  composition  of  wire  segment  N, 
for  N  =  1,2,. ..,Ns.  (Integer) 

=  1  implies  the  parameter  Loss  is  tan  Si  =  loss  tangent  of  wire  segment  N. 

=  2  implies  the  parameter  Loss  is  cr\  —  conductivity  of  wire  segment  N  (in 
0-1/  meter). 

=  3  implies  wire  segment  N  is  PEC,  and  DO  NOT  use  monopole  end-modes  on 
this  segment. 

=  4  implies  wire  segment  N  is  PEC,  and  DO  use  monopole  end-modes  on  this 
segment. 

Note  that  for  PEC  wire  segments,  the  user  should  only  set  Hoss(N)  =  4  when 
wire  segment  N  physically  connects  with  a  corresponding  wire  segment  in  an 
adjacent  lattice  cell.  This  enables  the  monopole  end-modes  to  account  for  con¬ 
tinuity  of  current  across  lattice  cell  boundaries  when  PEC  wire  segments  are 
used.  For  lossy  and/or  dielectric  wire  segments  (Hoaa(N)  =  1  or  2)  monopole 
end-modes  are  always  enabled. 

Radwire(N)  =  a  =  radius  (in  meters)  of  wire  segment  N,  for  N  =  l,2,...,Ns.  (Real) 

Epsrwire(N)  =  ejr  =  relative  permittivity  of  wire  segment  N, 
for  N  =  1,2,. ..,Ns.  (Real) 

Note  that  Epsrwire(N)  only  has  meaning  for  lossy  and/or  dielectric  wire  seg¬ 
ments,  i.e.,  when  Doss(N)  =  1  or  2. 

Loss  =  Loss  parameter  of  wire  segment  N,  for  N  =  1,2,. ..,Ns.  (Real) 

=  tan  Si  of  wire  segment  N  if  Hoss(N)  =  1. 
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=  <T\  of  wire  segment  N  (in  Sl~1/  meter)  if  Hoss(N)  =  2. 

Note  that  Loss  only  has  meaning  for  lossy  and/or  dielectric  wire  segments,  i.e., 
when  Hoss(N)  =  1  or  2. 

6.1.7  READ  7:  Lumped  Loads 

Hoad(N)  =  Location  of  lumped  load  N,  for  N  =  l,2,...,Nload.  (Integer) 

Zload(N)  =  Complex  impedance  of  lumped  load  N,  for  N  =  l,2,...,Nload. 
(Complex) 

Note  that  wire  location”  L  is  defined  as  follows: 

•  by  endpoint  A  of  segment  L  if:  L  <  Ns. 

•  by  endpoint  B  of  segment  (L  —  Ns)  if:  (Ns  +  1)  <  L  <  2  Ns. 

6.1.8  READs  8,  9  and  10:  Lattice  Geometry 

V(l)  =  v*  ss  r-component  of  the  lattice  defining  vector  v.  (Real) 

V(2)  =  vv  =  y-component  of  the  lattice  defining  vector  v.  (Real) 

W(l)  =  wK  sb  * -component  of  the  lattice  defining  vector  w.  (Real) 

W(2)  —  wy  =  y-component  of  the  lattice  defining  vector  w.  (Real) 

W(S)  =  w,  =  ^-component  of  the  lattice  defining  vector  w.  (Real) 

D(l)  =  Lattice  spacing  (in  meters)  in  the  u  =  i  direction.  (Real) 

D(2)  =  Lattice  spacing  (in  meters)  in  the  v  =  vx  x  + £  direction.  (Real) 

D(3)  =  Lattice  spacing  (in  meters)  in  the  w  =  v>x  x  +  f  +  to,  i  direction.  (Real) 

Note  that  the  lattice  defining  vectors  v  and  w  do  not  need  to  be  normalized  to 

unit  vectors  in  the  input  file,  only  the  direction  needs  to  be  defined.  Also  if  the  lattice 
is  perpendicular,  then  set: 
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•  V(l)  =  0, 

•  V(2)  =  1, 

•  W(l)  =  0, 

•  W(2)  =  0, 

•  W(3)  =  1. 

6.1.9  READ  11:  Direction  of  Propagation 

Uk(l)  =  x-component  of  propagation  direction  u*.  (Real) 

Uk(2)  =  y-component  of  propagation  direction  U*.  (Real) 

Uk(3)  =  x-component  of  propagation  direction  u*.  (Real) 

Note  that  the  propagation  direction  defining  vector  u*  does  not  need  to  be  nor¬ 
malized  to  a  unit  vector  in  the  input  file,  only  its  direction  needs  to  be  defined. 

0.2  Output  From  ADWIRS 

6.2.1  The  Output  File. 

The  program  ADWIRS  writes  its  output  file  to  standard  output  on  logical  unit  6. 
This  output  file  includes  such  information  as: 

1.  Input  data  and  any  errors  in  the  input  data. 

2.  Any  errors  encountered  during  the  computations. 

3.  The  impedance  matrix,  if  requested. 

4.  The  effective  wavenumber  he,  and  the  constitutive  parameters  as  calculated  by 
the  simple  formula  of  Equation  (6.1). 

5.  The  eigenfunction  currents,  dipole  moments  per  unit  cell  (P®,M#),  and  the 
average  electric  and  magnetic  fields  per  unit  cell  (E*,H°). 

6.  The  GPU  times. 
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6.2.2  The  Eigenfunction  Solution  File. 

The  program  ADW1RS  writes  the  eigenfunction  solutions  to  a  file  on  logical  unit 
9.  This  is  the  file  to  be  used  when  Equation  (6.1)  is  not  valid  and  the  polarization 
method  or  the  Maxwell’s  equations  method  is  to  be  implemented  by  the  user.  The 
data  can  be  read  by  the  following  FORTRAN  77  READ  statements. 

READ  (9,*)  ISVPTS , NUMMODES 
DO  I  -  1. ISVPTS 

READ  (9.*)  PARAM(I) ,0MEGA(I) 

DO  J  «  1 .NUMMODES 

READ  (9,*)  KE(I,J,1),KE(I,J,2) ,KE(I, J,3) 

READ  (9.*)  EPL(I,J,1) ,EPL(I,J,2) ,EPL(I,J,3) 

READ  (9.*)  HPL(I, J,l) ,HPL(I,J,2) ,HPL(I,J,3) 

READ  (9,*)  PO(I, J,l) ,PO(I, J,2) ,P0(I,J,3) 

READ  (9.*)  MO(I, J,l) ,MO(I, J,2) ,M0(I,J,3) 

END  DO 
END  DO 

The  above  variables  are  defined  here. 

NSWPTS  =  the  number  of  points  in  the  parameter  sweep. 

NUMMODES  =  the  number  of  MM  expansion  functions  on  the  wire  geometry. 

PARAM(I)  =  the  value  of  the  parameter  that  is  being  swept  at  sweep  point  I. 

OMEGA  (I)  =  the  angular  frequency  (in  radians  /  second)  at  sweep  point  I. 

KE(I,J,1)iKE(I,J,2),KE(I,J,S)  =  the  *,y,  z  components  of  the  effective  wave- 
vector  at  current  mode  solution  J  and  sweep  point  I. 

EPL(I,J,1),EPL(I,J,2),EPL(I,J»3)  =  the  x,y, z  components  of  E®  at  current 
mode  solution  J  and  sweep  point  I. 
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HPL(I,J,1),HPL(I,J,2),HPL(I,J,3)  =  the  x,y,z  components  of  H*  at  current 
mode  solution  J  and  sweep  point  I. 

P0(I,J,l),P0(I,J,2),P0(I,J,3)  =  the  x,y, z  components  of  P°  at  current  mode 
solution  J  and  sweep  point  I. 

M0(I,J,1),M0(I,J,2),M0(I,J,3)  =  the  x,y, z  components  of  M°  at  current  mode 
solution  J  and  sweep  point  I. 

6.2.3  The  Parameter  Sweep  Files 

When  a  parameter  sweep  is  performed  (either  a  standard  parameter  sweep,  or  a 
determinant  versus  he  sweep)  then  data  are  output  to  specific  files  for  plotting.  In  a 
standard  parameter  sweep  run,  tabular  listings  of  e^,  tan  6-  and  <r£  (as  computed  by 
Equation  (6.1))  versus  the  parameter  being  swept  are  written  to  files  on  logical  units 
41,  42  and  43,  respectively.  In  a  determinant  versus  kg  sweep,  tabular  listings  of  the 
absolute  value,  the  real  part,  and  the  imaginary  part  of  \Z\  versus  K/bo  are  written 
to  files  on  logical  units  20,  21  and  22,  respectively. 


96 


Chapter  7 
Summary 


An  artificial  medium  can  be  viewed  as  a  macroscopic  model  of  a  real  medium.  Typi¬ 
cally,  an  artificial  medium  consists  of  a  large  number  of  scattering  objects  distributed 
in  a  host  medium.  For  example,  the  artificial  media  treatable  by  the  solution  in  this 
dissertation  are  composed  of  a  3D  periodic  array  of  identical  arbitrarily-shaped  thin 
conductive  or  dielectric  wire  objects  arranged  in  the  homogeneous  host  medium.  In 
general,  the  wire  objects  perturb  the  eigenfunction  solution  for  a  plane  wave  in  the 
host  medium  such  that  a  plane  wave  with  a  different  effective  wavenumber  propagates 
in  the  artificial  medium.  Thus,  artificial  media  are  characterized  by  effective  consti¬ 
tutive  parameters,  and  in  general  are  anisotropic  media.  The  effective  constitutive 
parameters  can  be  a  function  of  frequency,  the  direction  of  propagation,  the  size, 
shape,  and  orientation  of  the  wire  objects,  and  the  constitutive  parameters  of  both 
the  wire  objects  and  the  host  medium. 

The  topic  of  this  dissertation  is  the  solution  of  the  plane  wave  that  propagates  in 
the  artificial  medium,  and  the  determination  of  the  effective  constitutive  parameters 
from  this  solution.  The  solution  is  obtained  by  formulating  an  integral  equation 
for  the  plane  wave  that  propagates  in  the  artificial  medium.  This  integral  equation 
is  solved  by  the  periodic  moment  method  (PMM),  and  yields  the  complex  effective 
wavenumber,  the  eigenfunction  currents  in  the  wire  objects,  and  the  eigenfunction 
fields  in  the  artificial  medium. 

Three  methods  are  presented  for  determining  the  effective  constitutive  parameters 
of  the  artificial  medium.  The  simplest  method  employs  the  effective  wavenumber  and 
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characteristic  impedance  of  the  artificial  medium,  and  applies  only  to  uniaxial  media 
with  propagation  and  polarization  along  principle  axes.  The  other  two  methods  apply 
to  general  anisotropic  media  at  any  direction  of  propagation.  One  method  enforces  the 
constitutive  relationship  equations  and  the  other  method  enforces  Maxwell's  source- 
free  equations  in  the  artificial  medium.  Both  methods  use  quantities  averaged  over 
the  volume  of  a  lattice  cell. 

The  solution  of  this  dissertation  has  several  distinguishing  characteristics  setting 
it  apart  from  other  solutions  to  artificial  media.  First  of  all,  no  static  approximations 
are  made,  and  the  only  limitation  on  frequency  is  that  the  wire  objects  and  lattice 
spacing  are  not  too  electrically  large.  This  solution  satisfies  Maxwell’s  equations,  in 
an  average  sense,  inside  the  artificial  medium.  Secondly,  this  solution  includes  the 
mutual  coupling  effects  of  the  3D  array.  Mutual  coupling  affects  the  fields  acting  on 
the  reference  wire  object  and  tl>e  current  shape  in  the  wire  objects.  Finally,  artificial 
media  composed  of  complex  wire  shapes  in  a  periodic  arrangement  can  be  analyzed. 

One  important  property  of  artificial  media  is  that  for  a  given  direction  of  propaga¬ 
tion,  there  are  two  distinct  plane  wave  modes  that  can  propagate  without  excitation. 
Note  that  for  some  artificial  media,  at  certain  directions  of  propagation,  one  or  both 
of  the  plane  wave  modes  may  be  the  same  as  a  plane  wave  in  the  host  medium. 
As  the  direction  of  propagation  changes,  the  plane  wave  modes  that  propagate  also 
change.  This  phenomenon  is  seen  by  observing  that  the  effective  wavenumber  for  a 
plane  wave  mode  of  propagation  is  a  strong  function  of  the  direction  of  propagation. 
This  is  also  true  for  real  anisotropic  media.  Also,  in  a  real  anisotropic  medium,  the 
permittivity  tensor  components  are  independent  of  the  direction  of  propagation.  It 
appears  that  the  effective  permittivity  tensor  components,  and  the  current  shape  on 
the  wire  objects,  can  be  a  function  of  the  direction  of  propagation.  Several  examples 
with  numerical  results  verify  these  observations. 
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Appendix  A 

The  Source  Pattern  Factor  for 
PWS  Current  Functions 


This  section  evaluates  the  source  pattern  factor  b)  given  by  Equation  (2.59) 

for  the  MM  expansion  functions.  The  axial  variation  along  the  length  of  the  straight 

■i 

current  elements  used  in  the  MM  expansion  functions  is  the  piecewise  sinusoidal  vari¬ 
ation  given  in  Equation  (2.18).  Inserting  this  current  variation  into  Equation  (2.59), 
and  integrating,  it  is  obtained  that 

= c~i±  -  <*(«)] 

where 

cr  =_ _ 1 _ 

nj±  ko  sin(fco  dnj)  [1  -  (a •  r±)2] 

A;i±(x)  =  •*■»***>  U(ini  •  f±)  «in(fco  *)  -  cos(*o  *)] . 

In  the  event  that  a  =  0  and  b  =  dnj,  as  in  the  Region  (I)  and  Region  (HI)  fields,  the 
source  pattern  factor  simplifies  to 

Kj±(*)  =  -1 

and 

K1±(b)  =  .  f±) rinffc, _  co>(fco 4^)1 . 


(A.1) 

(A.2) 

(A.3) 


A.l  Source  Pattern  Factor  Limits  of  Integration 


Figure  A.l  shows  the  three  primary  geometrical  arrangements  for  source  current  seg¬ 
ment  nj.  The  geometrical  arrangement  of  current  segment  nj  depends  upon  the  x 


Figure  A.l:  The  three  geometrical  arrangements  for  a  source  current  segment. 

values  of  its  endpoints,  which  dictate  the  expressions  for  the  electric  fields  and  the 
corresponding  limits  of  integration  for  the  source  pattern  factor,  P^j±(a,fe). 

The  arrangement  shown  in  Figure  A.l(a)  has  current  segment  nj  located  in  a 
plane  parallel  to  the  xy  plane.  Region  (11)  has  vanished  and  there  are  only  Region 
(I)  and  Region  (III)  type  fields.  Thus,  for  other  case,  the  limits  of  integration  for 
Pnj±(a,b)  are  a  =  0  and  6  =  d„j,  the  length  of  current  segment  nj. 

Similarly,  for  the  arrangements  shown  in  Figure  A.l(b)  and  (c),  if  the  field  point 
is  either  in  Region  (I)  or  Region  (Ill),  then  the  limits  of  integration  are  a  =  0  and 
b  =  dnj.  However,  if  the  field  point  is  in  Region  (II),  then  the  electric  and  magnetic 
fields  consist  of  both  right-going  and  left-going  waves,  as  given  in  Equations  (2.60) 
and  (2.63).  If  the  vector  from  the  origin  to  the  field  point  is  R,  then  the  point  V  on 
current  segment  nj  of  the  same  z  value  as  R  can  be  determined  by  solving 

i  •  R  =  £  •  (R,,,'  +  /' 

for  l',  resulting  in 

r=-^  +  —  z  =  a^+fi,it 

®nji  ®nj« 

where  R„jx  and  On,-,  are  the  x-components  of  R^-  and  respectively.  Employing 
the  notation  of  Equation  (2.60),  if  the  current  segment  geometry  is  as  shown  in  Figure 
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A.l(b),  then  the  limits  of  integration  are 


«nj+  =  0 

bnj+  =  anj  +  finj  S 
«nj-  =  Onj  +  fa  X 
bnj—  —  dn j. 

Similarly,  if  the  current  segment  is  as  shown  in  Figure  A. 1(c),  then  the  limits  of 
integration  are 

Onj-  =  0 

^nj-  =  ®n;  "I"  finj  31 

®nj+  =  ®nj  +  finj  * 

ftr »J+  =  dnj . 

A.2  Open-Current  End  Charge  Terms 

The  monopole  expansion  functions  have  a  discontinuous  current  distribution  in  which 
the  current  rises  sinusoidally  to  unity,  and  then  abruptly  falls  to  zero.  This  results  in  a 
charge  distribution  existing  at  the  open-current  end  of  monopole  expansion  functions. 
These  charge  terms  will  cancel  when  two  monopoles  are  placed  together  to  form  a 
dipole  expansion  function,  of  if  monopole  expansion  functions  are  used  to  model 
continuous  current  across  adjacent  lattice  cells.  In  both  these  cases,  the  physical 
current  is  continuous  and  does  not  contain  charges  resulting  from  an  open-current 
end. 

Since  the  geometry  of  interest  is  a  3D  periodic  array  of  wire  objects,  with  MM 
current  expansions  similar  (to  within  a  constant)  on  each  wire  object,  then  there 
exists  a  3D  array  of  charges  corresponding  to  the  open-current  monopole  expansion 
functions.  It  is  often  desirable  to  isolate  and/or  remove  the  contribution  to  the 
electric  field  from  this  3D  charge  array.  This  can  be  done  by  rearranging  the  product 
e„j±  P*i±( o,6),  which  occurs  in  the  expression  for  the  total  electric  field  of  Equation 
(2.69),  as 

=  c-it  {[B;i±(4)  +  Q*ni±(i)  -  BV±(.)  -  «:*(„)]  ?±- 
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where  and  A^j±(x)  are  given  by  Equations  (A.2)  and  (A.3),  respectively,  and 
Bnj±(*)  =  e^*0*^ •*±)  {-(&«, .  r ± )  cos(fco  *)  +  j  «in(*o  *)]  (A.5) 


and 


-]  «in(hp  X)  ikox(kmi  t4;) 

•  M  «  ■  V  C  • 


(A.«) 


q'  .  =  — i — 

w"’±l  '  *oC;J±  dn(»»<W 

The  term  (J'j±(*)  is  the  contribution  to  the  total  electric  field  from  the  charge  array 
located  at  endpoint  x.  Thus,  Qnj±(x)  can  he  omitted  to  remove  the  electric  field 
resulting  from  the  charge  array  at  endpoint  x.  Note  that  the  charge  contribution 
vanishes  when  *  =  0,  corresponding  to  zero  charge  at  the  zero  current  end  x  =  0. 
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Appendix  B 

Equivalent  Wire  Radius 


This  appendix  evaluates  ae,  the  equivalent  wire  radius,  as  defined  by  Newman  [25]. 
Basically,  the  mutual  impedance  between  two  filament  line  sources  separated  by  the 
equivalent  wire  radius  is  the  same  as  the  mutual  impedance  between  two  volumetric 
cylindrical  current  sources,  corresponding  to  the  thin  material  wires. 

Consider  a  2D  case  where  both  the  expansion  and  weighting  current  sources  are 
of  radius  a  and  are  oriented  along  the  x-axis,  with  S-directed  current.  The  radial 


variations  are  defined  as 


Je=CMk,p)  A  /  ma 


J«  =  — 

ir  a* 


A  /  ma 


where  the  constant  C  is  given  in  Equation  (2.17)  so  that  J,  and  Jw  have  unit  current. 
The  mutual  impedance  between  these  two  current  *ources  was  determined  by  Newman 
to  be 

-  ****>s8sl)  - 

(B.1) 

where  rfo  and  ho  are  the  characteristic  impedance  and  wavenumber  of  the  host  medium, 
and  may  be  complex.  Next,  Ztw  is  equated  to  the  mutual  impedance  between  two 
filament  line  sources  spaced  Oe  apart  from  each  other,  where  |hba«|  <  1,  i.e., 

Z,„  =  a.)]  (B.2) 
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This  expression  is  solved  for  o«,  resulting  in 


The  real  part  is  taken  since  o«  represents  a  real  physical  displacement. 


(B.3) 
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Appendix  C 

Stagger  of  Weighting  Monopole 


This  appendix  presents  the  method  of  staggering  the  weighting  monopole  segment  ms'. 
The  rules  for  staggering  the  weighting  monopole  segment  mi  depend  upon  the  rela¬ 
tive  geometrical  arrangement  of  both  monopoles  segments  mi  and  nj.  The  weighting 
monopole  segment  mi  is  staggered  by  an  equivalent  wire  radius,  a,,  in  a  direction  per¬ 
pendicular  to  its  centerline,  and  then  the  integration  is  performed  along  this  staggered 
monopole.  In  decreasing  order  of  priority  the  rules  are: 

1.  If  monopole  segments  mi  and  nj  are  both  in  the  same  plane,  and  that  plane 
is  parallel  to  the  xy  plane,  then  monopole  segment  ms'  is  staggered  a  distance 
of  tie  in  the  negative  x  direction.  This  stagger  is  not  absolutely  necessary  if 
monopole  segments  mi  and  nj  do  not  touch  or  overlap.  However,  this  stag¬ 
ger  will  increase  the  convergence  of  the  spectral  summation  due  to  the  alight 
distance  of  exponential  decay  for  the  evanescent  waves. 

2.  If  monopole  segments  ms  and  nj  are  collinear  and  touching,  and  are  parallel 
to  the  s-axis,  then  monopole  segment  ms'  is  staggered  a  distance  of  a*  in  the 
negative  x  direction. 

3.  If  monopole  segments  ms’  and  nj  are  collinear  and  touching,  but  are  neither 
parallel  to  the  s-axis,  nor  parallel  to  the  xy  plane,  then  monopole  segment  ms 
is  staggered  a  distance  of  Oe  in  the  unit  direction  of  the  vector  imi  x  (4«j  x  S). 
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4.  If  monopole  segments  mi  and  nj  are  touching,  bat  are  neither  colHnear,  nor  both 
parallel  to  the  xy  plane,  then  monopole  segment  mi  is  staggered  a  distance  of 
a,  in  the  nnit  direction  of  the  vector  i*,,  x  inj. 

5.  For  all  other  arrangements,  monopole  segment  mi  is  not  staggered  at  all,  and 
the  integration  is  performed  along  its  centerline. 
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Appendix  D 

Monopole  to  Monopole 
Arrangements 


As  shown  in  Figure  D.l,  there  are  eight  possible  geometrical  arrangements  for  the 
segment  mi  to  segment  nj  impedance  evaluation.  These  possible  arrangements  de¬ 
pend  upon  the  z  values  of  the  endpoints  of  the  segments,  denoted  as 
and  Znjff.  Conditions  on  the  endpoints  are  indicated  in  Figure  D.l,  where,  in  general, 
endpoint  b  has  a  greater  or  equal  z  value  than  endpoint  a.  The  additional  criteria  for 
choosing  which  case  the  segment  to  segment  impedance  evaluation  falls  under  can  be 
summarized  as: 


Case  1: 

Zmib  —  Zfija 

Case  2: 

^njb  ^  *mia 

Case  3: 

Zrnia  ^  %nja 

and 

Zjnib  ^  Znjb 

Case  4: 

*mia  ^ 

and 

*nju  ^  -  *nj6 

Case  5: 

*m»6  ^  *njb 

and 

Case  0: 

ia  ^  ^nja 

and 

?  *nifc 

Case  7: 

*mia  <  (*nja 

—  ^  *n »i6 

Case  8: 

*nja  ^  (*mia 
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Figure  D.l:  The  eight  potable  segment  mi  to  segment  nj  geometry  arrangements. 


Appendix  E 

The  Evaluation  of  Gmi  nj±(a,b) 


This  appendix  evaluates  Gmit„j±(a,  6),  defined  in  Equation  (3.32),  and  used  in  the 
evaluation  of  the  impedance  matrix  element  contribution  *  In  general,  Gmi>nj±(o,  b) 
requires  evaluation  when  at  least  a  portion  of  weighting  monopole  mi  is  within  Region 
(II)  of  expansion  monopole  nj.  To  repeat  Equation  (3.32),  Gmi,„j±(a,b)  is  defined  as 

=  [  (e„y±  • i».)  F^(l) 

The  function  P'i±  (a „j±,bnj±)  is  evaluated  in  Appendix  A.  Also,  the  scalar  product 
(e«»i±  •  ®m«)  is  incorporated  into  Gmiy„j±(a,  b)  so  that  the  contribution  from  open- 
current  end  charges  can  be  isolated  and/or  removed.  As  shown  in  Figure  E.l,  the 
limits  of  integration  (anj±,  6nj±)  depend  upon  the  position  vector  to  a  point  on  the 
weighting  monopole  mi,  and  hence  they  have  an  l  dependence.  For  a  given  value  of 
/,  the  integration  limits  are  determined  by  solving 

*  ‘  [Rni  +  I'injJ  =  *  •  [Rm,-  +  UU] 

for  resulting  in 

^  +  — »=«-.+ A-/  (e.i) 

®n;«  ®n;i 

where  Rmi,,  Rnjt,  «,»„  and  are  the  s-components  of  R*,,,  R^,  a*,,  and  a„„ 
respectively. 

As  was  shown  in  Figure  A.l(b)  and  A. 1(c),  there  axe  two  possible  arrangements 
of  monopole  nj  such  that  a  Region  (II)  field  exists.  This  results  in  two  different 
cases  of  integration  limits  for  the  njth  source  pattern  factor.  For  the  geometry  of 
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Figure  E.l:  Vector  relationship  for  the  center  Q  =  0  monopoles  mi  and  nj. 
Figure  A.l(b),  the  integration  limits  are 

®nj+  =  0 

bnj+  —  Ami  4"  Pmi  1 

®nj-  =  4*  Pmi  1 

bnj-  =  dnj 

and  for  the  geometry  shown  in  Figure  A.l(c)  the  limits  of  integration  are 


On j-  —  0 

bnj—  —  flmi  4"  Pmi  1 
®nj+  =  ®m<  4-  Pm  i  1 
bnj+  ~  ^nj  • 


110 


When  substituting  P'J±(anj±,^±)  into  the  defining  equation  for  Gm,tnj±(a,b),  the 
notation  of  Appendix  A  will  be  utilized,  as  follows: 


Figure  A. 1(b)  Geometry 

Figure  A. 1(c)  Geometry 

Pnj+(V>+) 

—  Pnj+(ami  +  Ami  0 

-®nj+(  V>+) 

~  ^*nj+ 

Bnj+(a*j+) 

=  C[ni+ 

■®ni+(ani+) 

=  Pnj+(a"»«  "b  Aw*  0 

Qnj+(Pnj+) 

—  Qnj+i^^mi  +  Ami  0 

Qnj+(.bnj+  ) 

=  Qnj+ 

Qnj+(anj+) 

=  Qnj+  =  0 

^nj+(anj+) 

~  Qnj+{ami  ”b  Ami  0 

A„j+(bnj+) 

=  Aj^.(Gtmi  +  Ami  /) 

■^n>+(^nj+) 

=  C2nj+ 

-^nj+(a»»J+) 

=  @2 nj+ 

Anj+  (anj+  ) 

—  ^nj+(®™i  "b  Ami  0 

Ki-fr nj-) 

=  C{nj. 

P»i-(6ni-) 

=  Bnj-  (otmi  +  Ami  0 

Bfij-  (anj—  ) 

=  Bnj-i01^*  "b  Ami  0 

Bnj-(anj~) 

=  Wnj- 

Qnj-{bnj~) 

=  Q“ 

'Vnj- 

Qnj-(bnj-) 

~  Qnj-[afni  "b  Ami  0 

Qnj-(anj-) 

=  Qnj-  (°™i  "b  Ami  0 

Qnj-  ( anj -  ) 

=  Q'nj-  =  0 

Anj-  finj-  ) 

=  Cfc*. 

Kj-iKj-) 

=  A*j_(o^„,-  +  Ami  0 

A„j-  (anj-  ) 

=  A*J_(amj  +  Ami  1) 

Anj-(anj~) 

=  <?*„>-• 

(E 

The  terms  C'nj± 

and  Cjnj±  are  constants  with  respect  to 

z  which  result  from  t 

endpoints  of  segment  nj.  Similarly,  the  terms  Q„}±  are  the  charge  terms  associated 
with  the  open-current  at  the  endpoints  of  segment  nj,  and  are  independent  of  /.  The 
charge  terms  (<**».  +  Am*  /)  must  be  removed  from  the  integral  for  the  answer  to 
be  right,  but  I’m  not  sure  why  this  is  so.  Substituting  Fmi(l)  and  the  above  notation 
for  Pnj±(anj±,bnj±)  into  the  defining  equation  for  Gmiitlj±(a,b),  and  rearranging, 

Gmi,nj±  (®J  b)  =  S  Cnj  [  C\±  I\±  +  Ci±  ^2±  +  c>3±  ^3±  ] 

*** 

where  S  is  given  in  Equation  (4.19),  and 

C'±  =  ((*"»• '  *±)  “  i(®"»  •  *»i) (*"J  *  *±)1 

°2±  =  sinfodni)  “  *(*"“'  '  f±)  ‘  f±)l 

=  Cjnj±  (firo,-  •  &nj)  —  {p\nj±  +  Qnj±)  (**»•  *  *±) 
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K±  —  eitoamifa, jt±) 

h±  =  ^*exp  -  P™*nj]  •  f±)  sin  [fco(ot»<  +  /3mi  l)]  si n(M)  dl 

h±  =  /  exp  (-iMa*,.  -  &•»•»,]  •  r±) cos  [fco(ami  +  0mi  l)J  nn(*o/)  dl 

*a 

h±  =  J  exp  [—jkol^Bmi  •  r±)]  sin(fcol)  dl. 

Note  that  the  definition  of  £7a±  contains  the  open-current  end  charge  term  Qnj±- 
This  term  can  be  omitted  when  two  monopoles  are  placed  together  to  form  a  dipole 
expansion  function.  In  this  case  the  charges  from  the  two  monopoles  will  cancel 
exactly. 

In  the  evaluation  of  I\±  and  h±  it  is  worthwhile  to  define 

^niidt  =  jk)[*mi  Ani®nj]  * 

/ 

Also,  the  trigonometry  identities 


sin  (A  +  B)  =  sin  A  cos  B  +  cos  A  sin  B 

cos(  A  +  B)  —  cos  A  cos  B  —  sin  A  sin  B 
2  sin  A  sin  B  =  —  cos  (A  +  B)  +  cos(  A  —  B) 
2  cos  A  sin  B  —  sin  (A  +  B)  —  sin(  A  —  B) 
are  applied  to  the  Ii±  and  Ii±  integrands  to  obtain 


h±  =  i  [/i+  -  /ri  -  5  in*  -  /n 

hi  =  i  cosffcoow)  [JJt  -/f)  +  i  C08(fcoam.)  [/Jf  -  /J- J 


where 

n* 


-jf' 


«in[MA».  +  l)l\  dl  = 


(«m,±)2  +  (M/9U.  +  l)]2 

X  (~wmi±  sin[*o (0mi  +  1)6]  -  *o(/9m,  + 1) cos[*o(/Sm,  +  1)6])  - 

-  e-Wmi±a  (-wm,±  sin[ho(/9m,  +  l)o]  -  k>(A»*  +  1)  cos[*o(/0m,  +  l)a])] 

1 


-i 


swlh(0m,-l)ll<U  = 
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x  [e—**  (-uw±  «m[*o(A».  -  1)4]  -  k,(0„.  -  1)  cm[*o(A..  -  1)4])  - 
-  (-wmi±  rin[*o(A~  -  1W  -  ko(A«i  -  1)  co»(fe>(ft»j  -  l)a))] 

*  =  -w*- + IM<“  = 

X  (ko(P„ j  +  1)  >m(4o(^.  +  1)4]  -  UWt  CO.fc(A.j  +  1)4])  - 

-  (*o(flni  +  1)  >in[*o(A»,  +  1)«]  -  “mi±  CO«[*o(A»i  +  1)<*1)] 

/r  =  -1),|<“  =  &+1wlTI)Fx 

x  (*o(Am  -  1)  «in(fc(A*  -  1)4]  -  wmj±  co.MA.i  -  1)4])  - 

-  (ko(/3„i  -  l)«m[4o(/3m,  -  l)o]  -  iK„i±  co»[4o(A»i  -  l)a])] 

The  above  forms  of  1\±  and  1?±  axe  valid  as  long  as  /9mi  /  ±1,  in  which  case 


H*  =  ii'  =  o 


always 


_JC~_  b~a  tf  wm.±  =  I 

~  *  ~  otherwise. 

— Wmii 


The  J3±  integral  can  be  evaluated  directly  as 


h±  =  Ja  exp  ■  f±)]  sin(M)  dl  =  ^  ?±)2] 

X  [e-iMCW*)  (_i(ifji. .  f±)  Bin(fco6)  _  c08(Ao6))  - 
-  e_jtoa(*m,"f±)  •  r±) sin(fcoa)  -  cos(fea))  -] . 
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Appendix  F 

Integration  Limits  for  Evaluating 

Zz=:  ■ 

mi,nj 


This  appendix  defines  the  limits  of  integration  over  the  expansion  monopole  nj  and 
weighting  monopole  mi,  as  needed  in  the  evaluation  of  the  impedance  term 
In  general,  this  includes  the  integration  limits  for  the  following  terms: 

Pnj-  (°»i-  »  bnj-  )  (°m»l  >  ^mil )» 

^nj+  (®nj+ 1  (®n>«3>  bmo) 

(*mi,nj±  (®m»2  j  hmia)* 

Recall  from  Figure  D.l  that  there  are  eight  possible  geometrical  arrangements  for  the 
monopole  mi  to  monopole  nj  mutual  impedance  contribution  Z^nj.  Furthermore, 
as  shown  in  Figure  A.l,  there  are  three  primary  arrangements  for  each  of  monopoles 
mi  and  nj.  Each  case  is  analyzed  in  greater  detail  below. 


Case  1: 


As  shown  in  Figure  D.l  for  a  Case  1  geometry,  monopole  mi  is  entirely  in  Region  (I) 
of  monopole  nj.  The  integration  limits  are 


Integration  limits: 


®rnil  —  0,  &m*l  —  4ni 
|  ®nj-  =  bnj-  —  &nj 


and  then  (««,-, Kj-)  is  evaluated  as  required  in  Equation  (3.26). 
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Case  2: 


As  shown  in  Figure  D.l  for  a  Case  2  geometry,  monopole  mi  is  entirely  in  Region 
(III)  of  monopole  nj.  The  integration  limits  are 


Integration  limits: 


®mi3  ~  0,  imi3  =  dm. 

< 

Onj-t-  =  0,  bnj+  =  <Lnj 


and  then  Pnj+  (**»*+ » bnj+  )  Pmi+  (**»»! »  fcn.ii)  is  evaluated  as  required  in  Equation  (3.27). 


Case  3? 

As  shown  in  Figure  D.l  for  a  Case  3  geometry,  monopole  mi  is  entirely  in  Region  (II) 
of  monopole  nj.  The  integration  limits  are 

Integration  limits:  j  ami3  =  0,  fcml2  =  dm. 

and  then  Gmittlj±(ami2,  fcmi2)  is  evaluated  as  required  in  Equation  (3.31). 


Case  4: 


As  shown  in  Figure  D.l  for  a  Case  4  geometry,  monopole  mi  is  within  both  Region 
(I)  and  Region  (II)  of  monopole  nj.  Figure  F.l  shows  the  four  different  possible  ways 
that  the  Case  4  geometrical  arrangement  can  occur.  For  the  integration  limits  that 
follow,  it  is  necessary  to  define  the  following  terms: 


(F.l) 


Alim  _ 

"•»'  cimi. 

D lim  _  RnJi~Rmi*  i  anii  J  _ 

•D"»*  “  ami.  +  a„i,°"J- 

The  integration  limits  for  the  evaluation  of  P*^_(any_,fcnj_),  Pm,_(amii,fcmi1)  and 
£mi,ni±(om>2*6ii»2),  u  required  by  Equation  (3.33),  are  given  in  Table  F.l. 


Case  5: 

As  shown  in  Figure  D.l  for  a  Case  5  geometry,  monopole  mi  is  within  both  Region 
(III)  and  Region  (II)  of  monopole  nj.  Figure  F.2  shows  the  four  different  possible  ways 
that  the  Case  5  geometrical  arrangement  can  occur.  The  integration  limits  for  the 
evaluation  of  P*j+(anj+,bnj+),  P;,l+ («„,*,&*,*)  and  Gml>ji:(ami2,fcmi2),  as  required 
by  Equation  (3.33),  are  given  in  Table  F.2. 
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Case  6: 


As  shown  in  Figure  D.l  for  a  Case  6  geometry,  monopole  mi  is  within  all  Regions 
(I), (II)  and  (III)  of  monopole  nj.  Figure  F.3  shows  the  four  different  possible  ways 
that  the  Case  6  geometrical  arrangement  can  occur.  The  integration  limits  for  the 
evaluation  of  P'y_ (o^_ ,  6^- ),  P^,_  (omit ,  &m.i),  P„*j+  (<*„;+  ,&«;+),  P^i+  («™a,  &mis) 
Gmt,nj±(omrt, bmrt))  u  required  by  Equation  (3.35),  are  given  in  Table  F.3. 


Case  7; 

As  shown  in  Figure  D.l  for  a  Case  7  geometry,  Region  (II)  of  monopole  nj  has 
vanished,  and  monopole  mi  is  within  Region  (I)  and  Region  (III)  of  monopole  nj. 
Figure  F.4  shows  the  two  different  possible  ways  that  the  Case  7  geometrical  ar¬ 
rangement  can  occur.  The  integration  limits  for  the  evaluation  of  P„j_(anj-,bnj-), 
Pmi-(°mii,&mii)»  *%+(*«;+,  &*,+)  i£,-+(«mi3.*b«s),  “  required  by  Equation  (3.36), 

are  given  in  Table  F.4. 

Case  8; 

As  shown  in  Figure  D.l  for  a  Case  8  geometry,  monopole  mi  is  entirely  within  Region 
(II)  of  monopole  nj.  However,  monopole  mi  has  no  extent  in  the  z  direction,  and 
Cfniliy±(am.2»5m.2)  decouples,  as  given  in  Equation  (3.37).  Figure  F.5  shows  the  two 
different  possible  ways  that  the  Case  8  geometrical  arrangement  can  occur.  The 
following  term  is  required  for  presenting  the  integration  limits: 

,  Pf nit  ■&! j* 

“njO  =  - -  ’ 

fliy* 

The  integration  limits  for  the  evaluation  of  Pnj±(anj±ibnj±)  and  ^,±(«m<2,^ *«)»  u 
required  by  Equation  (3.36),  are  given  in  Table  F.5. 
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J 

Figure  F.l:  The  four  different  Case  4  impedance  possibilities. 


Table  F.l:  Integration  limits  for  a  Case  4  impedance. 


'/ 

figure  F.2:  The  four  different  Case  5  impedance  possibilities. 


Table  F.2:  Integration  limits  for  a  Case  5  impedance. 
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Figure  F.3:  The  lour  different  Cue  6  impedance  possibilities. 


Table  F.3:  Integration  limits  for  a  Case  6  impedance. 
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Figure  F.4:  The  two  different  Case  7  impedance  possibilities. 


Table  F.4:  Integration  limits  for  a  Case  7  impedance. 
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Figure  F.5:  The  two  different  Case  8  impedance  possibilities. 
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Table  F.5:  Integration  limits  for  a  Case  8  impedance. 
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Appendix  G 

The  Ellipsoid  of  Wave  Normals 


Consider  a  homogeneous  real  anisotropic  non-magnetic  medium  arranged  such  that 
the  coordinate  axes  are  aligned  with  the  three  axial  directions  (called  the  principle 
dielectric  axes)  x,y  and  x.  It  is  shown  by  Sommerfeld  [23,  Sec.  4.24]  that  any  real 
anisotropic  medium  can  be  arranged  in  such  a  way.  This  is  a  uniaxial  medium  with 
permittivity  tensor  given  by 


«e 


*,0  0 

0  e,  0 

0  0  c, 


where  e„  eB  and  e,  are  called  the  principle  dielectric  constants.  The  fields  D  and  E 
are  related  as 


D,  =  CgEg 


Dy  —  CyEy 

Dt  =  e„Et 


The  electric  energy  density  is  given  by 

1  _  _  1  (Dt  Dl  I%\ 

«,.  =  -E-D  =  -(— +  -J  +  — J. 


Letting  C  =  8tu>„  and  writing  z,y  and  x  in  place  of  DK/y/CtDy/y/C  and  DM/y/C, 
and  considering  these  as  Cartesian  coordinates  in  space,  we  get  the  equation  of  the 
e&pxoid  of  wave  normals 


x*  y2  x* 

—  +  —  +  —  =  1. 
6,  ty  e. 


(G.1) 
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This  ellipsoid  can  be  used  to  predict  the  two  phase  velocities,  and  hence  the  two 
wavenumbers  of  a  plane  wave  propagating  in  any  direction  u*,  as  follows  [24,  Sec. 
14.2.3].  Draw  a  plane  normal  to  u*  through  the  origin.  The  intersection  of  this 
plane  with  the  ellipsoid  is  an  ellipse.  The  major  and  minor  semi-axes  of  this  ellipse 
are  proportional  to  l/vp,  the  reciprocal  of  the  phase  velocities.  Once  the  two  phase 
velocities  have  been  determined,  the  corresponding  wavenumbers  can  be  found  as 


where  u  is  the  angular  frequency  in  radians  per  second. 

These  results  can  be  applied  to  artificial  media  in  the  following  manner.  Consider 
an  anisotropic  uniaxial  artificial  dielectric.  Assume  that  the  effective  wavenumber 
corresponding  to  one  of  the  two  roots,  has  been  determined  for  propagation  along 
two  of  the  principle  dielectric  axes.  For  example,  for  one  of  the  two  roots,  assume 
he  =  h,  for  propagation  along  the  x-axis  and  he  =  ky  for  propagation  along  the  y- 
axis.  We  wish  to  compute  the  effective  wavenumber  for  propagation  in  any  direction 
in  the  xy-plane,  i.e.,  u*  =  xu*,  +  y  u*y. 

See  Figure  G.l  for  the  geometrical  representation  of  the  solution.  The  ellipse 
equation  can  be  written  as 

x 2  y 2 

+  =  (G.2) 


Now,  substituting  x  —  a  and  y  =  b  =  a  tan  a  =  a(titv/u*x)  into  Equation  (G.2),  and 
solving,  it  is  obtained  that  for  the  given  direction  of  propagation  at  angle  a, 


he  =  Va2  +  b?  =  kx  ky 


y/l  +  tan2  a 
/ft2  -f  k2  tan3  a 


(G.3) 
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Appendix  H 
Material  Spheres 


This  appendix  presents  an  integral  equation  and  periodic  method  of  moments  (PMM) 
[2,  3,  4]  solution  to  the  problem  of  determining  the  effective  permittivity  of  an  arti¬ 
ficial  dielectric  composed  of  a  3D  periodic  array  of  small  identical  dielectric  material 
spheres.  Note  that  the  material  presented  in  this  appendix  is  mostly  self-contained, 
i.e.,  there  are  few  references  to  material  in  this  dissertation  outside  this  appendix. 
However,  much  of  the  theory  is  similar  to  that  presented  in  Chapter  2,  and  thus 
there  is  some  repetition.  The  major  differences  are  those  necessary  to  account  for  the 
scattering  objects  being  material  spheres  rather  than  wire  objects.  Also,  field  com¬ 
putations  are  made  by  employing  approximations  in  the  space  domain,  as  opposed  to 
the  spectral  domain  methods  presented  in  Sections  2.4  and  2.5. 

As  with  the  method  presented  in  Chapter  2,  this  method  is  based  upon  finding  the 
complex  wavenumber  for  a  plane  wave  propagating  in  the  artificial  dielectric,  from 
which  the  complex  effective  permittivity  can  be  deduced.  As  before,  mutual  coupling 
between  the  small  spheres  is  included  in  the  PMM  formulation. 

The  scattering  objects  are  identical  dielectric  spheres  arranged  in  a  3D  periodic 
lattice.  When  a  plane  wave  propagates  through  an  artificial  dielectric,  currents  are 
induced  on  or  in  the  scattering  objects.  These  currents  can  be  viewed  as  macro¬ 
scopic  current  moments,  analogous  to  the  microscopic  dipole  moments  induced  in  the 
molecules  of  an  actual  dielectric  [1].  The  effect  of  the  macroscopic  current  moments 
is  to  produce  a  net  current  moment  per  unit  volume,  and  thus  the  artificial  dielectric 
has  some  complex  effective  permittivity. 
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In  the  present  case,  the  complex  effective  permittivity  of  the  artificial  dielectric 
is  a  function  of  the  electrical  size  and  spacing  of  the  dielectric  spheres,  and  the  ma¬ 
terial  parameters  of  both  the  host  medium  and  the  spheres.  By  properly  choosing 
these  values,  it  may  be  possible  to  design  an  artificial  dielectric  medium  of  desired 
permittivity  and  loss  tangent.  A  recent  application  of  artificial  dielectrics  is  in  the 
microwave  welding  of  polymers  [11].  In  this  case,  a  lossy  dielectric  of  desired  con¬ 
ductivity  is  produced  by  the  proper  mixture  of  HC1  doped  polyaniline  particles  in  a 
polyethylene  host. 

Section  H.l  of  this  paper  presents  the  integral  equation  of  the  artificial  dielec¬ 
tric  medium.  This  integral  equation  is  solved  by  the  periodic  moment  method  using 
expansion  functions  suitable  for  small  dielectric  spheres.  Details  concerning  the  eval¬ 
uating  of  the  impedance  matrix  are  given  in  Section  H.2.  Numerical  results  are 
presented  in  Section  H.3,  illustrating  the  complex  effective  permittivity  for  different 
artificial  dielectric  compositions. 

H.l  Theory 

H.1.1  Derivation  of  the  Integral  Equation 

This  section  will  present  the  integral  equation  and  PMM  solution  for  a  plane  wave 
propagating  through  an  artificial  dielectric  medium  composed  of  small  dielectric 
spheres  arranged  on  a  periodic  lattice.  This  solution  will  yield  the  complex  wavenum¬ 
ber  of  the  plane  wave,  and  in  turn,  the  complex  effective  permittivity  of  the  artificial 
dielectric.  Throughout  this  appendix,  all  fields  and  currents  are  assumed  to  be  time 
harmonic  with  the  e7"*  time  dependence  suppressed. 

As  shown  in  Figure  H.l,  the  geometry  of  the  artificial  dielectric  consists  of  a 
3D  triple  periodic  array  of  identical  dielectric  spheres  located  in  a  homogeneous  and 
isotropic  host  medium.  The  spheres  have  radius  a,  and  are  homogeneous  with  per¬ 
meability  and  permittivity  denoted  by  (/io,Ci),  and  wavenumber  kx  =  Wy/poeT.  The 
homogeneous  host  medium  has  permeability  and  permittivity  denoted  by  (po«€o)i 
wavelength  Ao,  and  wavenumber  ko  =  Uy/fHo^o.  Note  that  this  host  medium  is  not 
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Figure  H.l:  Geometry  of  a' 3D  periodic  artificial  dielectric  composed  of  dielectric 
spheres. 

necessarily  free  space  and  may  be  lossy.  The  spheres  are  arranged  in  a  rectangular 
lattice  with  spadngs  ds,  dy  and  d„  in  the  i,  y  and  i  directions,  respectively.  They  are 
referenced  by  the  index  Q  =  (ix,jv1kM)  where  —  oo  <  (*x,jy,  Jfe,)  <  oo.  The  reference 
or  center  sphere  is  centered  at  the  origin  and  is  indexed  by  Q  =  0  =  (0,0,0).  Also, 
let  AQ  =  ixdxx  + jvdyy  +  kxdx&  be  the  position  vector  from  the  origin  to  the  center  of 
sphere  Q.  Typically  there  are  a  large  number  of  electrically  small  spheres  per  cubic 
wavelength. 

Following  the  general  methods  presented  by  Blanchard  and  Newman  [21, 22],  the 
effective  permittivity  of  the  artificial  dielectric  medium  is  determined  by  first  ■■■n*ning 
that  a  plane  wave  of  unknown  wavenumber  is  propagating  through  the  medium.  If 
the  plane  wave  is  propagating  in  the  u*  direction,  then  it  will  have  spatial  variation 
of  the  form 

«-*•*,  (H.l) 

where  ke  =  he  U*  is  the  unknown  complex  wave-vector,  and  R  =  xk  +  yy  +  ri  is  a 
position  vector  from  the  origin  to  the  point  (x,y,z).  For  a  given  direction  u&,  it  is 
desired  to  find  he  such  that  this  plane  wave  satisfies  Maxwell's  source  free  equations 
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and  all  of  the  boundary  conditions  in  the  artificial  medium,  i.e.,  the  normal  mode  of 
propagation  for  the  artificial  dielectric  medium.  Once  the  propagation  constant  fc* 
is  known,  then  the  effective  permittivity  of  the  artificial  dielectric  medium  is  found 
through  the  relationship 

k3 

K  =  Wy/iOte  or  e,  as  — (H.2) 

where  e*  is  the  complex  effective  permittivity  of  the  artificial  dielectric  medium.  Note 
that,  in  general,  the  artificial  dielectric  is  an  anisotropic  medium,  with  different  di¬ 
rections  Ufc  yielding  different  values  for  K  and  e*.  However,  since  the  spheres  are 
symmetric  and  closely  spaced,  the  directional  variation  of  e*  should  be  negligible. 

In  formulating  the  integral  equation  for  the  artificial  dielectric  medium,  the  volume 
equivalence  theorem  is  used  to  replace  the  dielectric  spheres  by  the  host  medium  and 
the  equivalent  electric  volume  polarization  currents  [1,  Sec.  7.7] 

J  =  M*i  -CojE3,  (H.3) 

where  E1  is  the  total  electric  field  inside  the  dielectric  spheres.  As  illustrated  in 
Figure  H.2,  the  dielectric  spheres  are  replaced  by  the  host  medium  and  the  equivalent 
volume  currents  J.  This  current  J,  which  exists  in  each  and  every  sphere,  is  written 
as 

J  =  £  (H.4) 

Q 

where  J**  is  the  current  in  the  volume  occupied  by  sphere  Q.  In  Equation  (H.4) 
and  others  to  follow,  it  is  implicit  that  the  summation  is  over  all  Q,  i.e.,  — oo  < 
(*x,ju,kj)  <  oo.  Since  we  seek  a  solution  to  Maxwell’s  source  free  equations,  there 
are  no  impressed  currents,  and  thus  E*  is  the  electric  field  of  J  radiating  in  the 
homogeneous  host  medium.  Equation  (H.3)  can  be  rearranged  as  the  homogeneous 
equation 

-  E*  +  t~7~~ - r  as  0  in  each  V**,  (H.5) 

-  e<>) 

and  is  to  be  solved  for  the  complex  wavenumber  ke,  and  the  current  in  the  center 
Q  =  0  sphere,  by  the  PMM. 
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Figure  H.2:  Equivalent  representation  for  the  3D  artificial  dielectric  of  spheres. 


Due  to  the  periodic  nature  of  the  array  of  spheres,  and  of  the  plane  wave  of 
Equation  (H.l),  the  current  is  identical  in  each  sphere  except  for  an  amplitude  and 
phase  change  corresponding  to  the  amplitude  and  phase  of  the  plane  wave  at  the 
center  of  the  sphere.  In  other  words,  the  current  in  sphere  Q  differs  from  the  current 
in  the  center  or  reference  sphere  by  the  complex  multiplier 

CQ  =  e-A  AQ  (H.6) 


As  a  result,  the  only  unknowns  are  K  and  the  current  in  the  center  or  reference 
sphere. 

H.l. 2  PMM  Solution  of  the  Integral  Equation 

The  first  step  in  the  PMM  solution  is  to  expand  the  unknown  current  J  as 

J  =  (H.7) 

Q  Q  n=l 

where  the  J?  are  N  linearly  independent  expansion  functions  for  the  current  in  sphere 
Q,  and  the  1„  are  N  unknown  expansion  coefficients,  with  n  =  1,2, . . .  ,N.  Due  to 
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the  periodic  nature  of  the  problem,  it  is  only  necessary  to  match  Equation  (H.5)  over 
V9,  the  volume  of  the  center  sphere.  Next  define  N  linearly  independent  weighting 
functions  in  the  center  sphere,  denoted  as  Wm  with  m  =  1,2, ...,  N.  Substituting  J  of 
Equation  (H.7)  into  Equation  (H.S),  and  taking  the  inner  product  of  the  result  with 
the  N  weighting  functions,  reduces  Equation  (H.5)  to  an  order  N  matrix  equation 
which  can  be  written  as 

[Z  +  AZ]  1  =  0.  (H.8) 

Here  [Z  +  A Z]  is  the  order  N  impedance  matrix  and  I  is  the  length  N  solution  vector 
containing  the  /„  expansion  coefficients  of  Equation  (H.7).  The  impedance  matrix 
elements  are  given  by  (m,n  =  1,2,...,  JV) 

z„„(K)  =  £  c<*zS„  =  -  £  <?Q  /  E?  ■  Wm  iv  (H.9) 

Q  Q  m 

1 

where  is  the  electric  field  of  Jj3,  the  nth  expansion  function  in  sphere  Q,  radiating 
in  the  homogeneous  host  medium.  The  integration  limits  Vm  correspond  to  the  region 
occupied  by  Wm.  Although  the  A  Zmn  do  not  depend  on  K,  the  Zmn  do,  because  of 
their  dependence  upon  C^.  Note  that  the  A Zm„  terms  are  non-zero  only  when  the 
expansion  functions  Jj^  are  in  the  reference  Q  =  0  sphere.  Section  H.2  discusses  the 
numerical  evaluation  of  the  impedances  in  Equations  (H.9)  and  (H.10). 

The  homogeneous  matrix  Equation  (H.8)  will  have  a  non-trivial  solution  only  if 
the  determinant  of  the  impedance  matrix  is  zero.  Thus,  he  is  found  by  solving  the 
fundamental  equation 

|Z(*e)  +  AZ|  =  0,  (H.ll) 

usually  on  an  iterative  basis. 

H.1.3  The  MM  Expansion  and  Weighting  Functions  for 
Small  Spheres 

In  this  section,  MM  expansion  and  weighting  functions  suitable  for  electrically  small 
spheres  are  defined.  For  spheres  of  arbitrary  size,  the  sphere  eigenfunctions  [37] 
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form  a  reasonable  set  of  expansion  functions.  However,  for  electrically  small  spheres 
(a  <  Ao),  the  shape  of  the  current  in  a  given  sphere  is  very  close  to  the  lowest  order 
eigenfunction  current  for  an  isolated  sphere  of  the  same  radius  and  permittivity  in 
the  homogeneous  host  medium.  Essentially,  the  higher  order  terms,  which  drop  off  as 
higher  powers  of  |Aoa|  are  being  ignored,  and  only  the  lowest  order  eigenfunction  cur¬ 
rent  is  retained  in  each  sphere.  A  Galerkin  PMM  solution  is  employed  with  weighting 
functions  chosen  identical  to  the  expansion  functions  in  the  center  sphere.  For  N  =  1 
expansion  function,  the  fundamental  equation,  Equation  (H.ll)  reduces  to 

Zn(*e)  +  AZn  =  0,  (H.12) 


i.e.,  the  self  impedance  of  the  single  expansion  function  is  zero. 

The  expansion  function  is  defined  by  illuminating  the  isolated  center  sphere  with 
an  x  polarized  plane  wave  propagating  in  the  +z  direction  (u*  =  z),  and  with  the 
incident  electric  field 

E1  =  xe~ikoM.  (H.13) 

J® ,  the  n  =  1  expansion  current  in  the  Q  =  0  sphere,  is  the  current  induced  by  this 
incident  field  in  the  isolated  center  sphere.  This  current  is  given  by  [37] 

Jf  =  C,  [aim*!',  +  iiing',] ,  (H.14) 

where  a\  and  6!  are  the  eigenfunction  expansion  coefficients,  and  mg'  and  ng,  are 
the  appropriate  spherical  vector  wavefunctions  inside  the  sphere.  C\  is  an  arbitrary 
normalization  constant  chosen  here  such  that 


Jy0  dv  =  Cif  §  [ojmili  +  jb[n£\]  dv  =  lx  (Am). 

The  eigenfunction  coefficients  can  be  written  as 

at  _  ji(po)hj2)(po)  ~  M2)(Po)jo(Po) 

bt  = _ M2)(po)io(Po)  ~  ii(Po)^o2)(Po) _ 

1  M2)(po)ji(pi)  -  «ii(pi)42)(po)  +  (l*=l)  ii(pi)M2)(po)’ 


(H.15) 


(H.16) 

(H.17) 
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where  jm(  )  end  A^(-)  ere  n**  order  iphericel  Bessel  end  sphericel  Henkel  functions 
of  the  second  type,  respectively,  po  =  koa,  p\  =  hja  end  *  =  ei/co.  Denoting  p  =  k\t 
where  r  is  the  rediel  distance  from  the  center  of  the  sphere,  the  sphericel  vector 
wevefunctions  inside  the  sphere  can  be  written  es 

m„ii  =  ji(p)  [cos ^6  —  cot 6 nn 4 (H.18) 

^  [2ji(p)  sin  9  cos  ^  #  +  (pji(p)]'  cos  9  cos  4  $  +  (pji(p)]'  sin  4  #] ,  (H.19) 

where  the  prime '  denotes  differentiation  with  respect  to  p  =  fcir,  end  6  end  ^  ere  the 
usual  spherical  angle  coordinates.  Note  that  the  expansion  function  J?  in  sphere  Q  is 
identical  to  the  expansion  function  J*  in  the  center  sphere  except  for  the  translation 
AQ. 


H.2  Evaluation  of  the  Impedance  Matrix 

This  section  describes  the  numerical  evaluation  of  the  impedance  matrix  element 
[Zu  +  A Zn]  of  Equations  (H.9)  and  (H.10)  using  the  expansion  and  weighting  func¬ 
tion  given  in  Section  H.1.3.  For  simplicity,  the  direction  of  propagation,  fi*  of  Equa¬ 
tion  (H.l),  is  assumed  to  be  in  the  yx  plane.  Thus,  the  complex  multiplier  of  Equa¬ 
tion  (H.6)  does  not  depend  on  the  x  lattice  spacing  index  »„  and  is  written  as 

qQ  —  jwd9<x»++k.d,un1>)'  (H.20) 

where  i>  is  the  angle  between  ft*  and  the  y-axis. 

As  indicated  by  Equation  (H.9),  Zu  is  obtained  by  summing  Z$,  weighted  by  the 
multiplier  C**,  for  every  expansion  current,  J?,  in  the  infinite  3D  periodic  array.  Z^ 
is  the  mutual  impedance  between  expansion  function  J?  in  sphere  Q  and  weighting 
function  Wi  =  J*  in  the  center  sphere.  In  order  for  the  triple  sum  of  Equation 
(H.9)  to  converge  absolutely,  the  mutual  impedances,  Zff,  must  drop  off  with  radial 
distance  r  more  rapidly  than  1/r3.  This  is  not  the  case,  since  for  large  r  the  field  E? 
(and  thus  Z$)  drops  off  only  as  1/r.  However,  due  to  cancellations  caused  by  the 
oscillatory  behavior  of  C**  and  the  triple  sum  does  converge,  although  its  rate 
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of  convergence  is  very  slow.  Thus,  some  method  is  desired  to  quickly  evaluate  this 
slowly  convergent  summation. 

Many  techniques,  based  upon  mathematical  identities  such  as  the  Poisson  sum 
formula,  have  been  presented  to  accelerate  the  convergence  of  a  series  [20,  38].  Here 
we  present  a  method  which  is  based  more  upon  physical  reasoning  than  mathematical 
identities.  The  method  is  based  upon  the  assumption  that  the  dielectric  spheres  are 
electrically  very  small  and  are  on  a  lattice  spacing  which  is  electrically  very  small.  The 
method  employs  increasingly  crude  approximations  as  the  distance  from  sphere  Q  to 
the  center  sphere  increases.  As  illustrated  in  Figure  H.3,  the  spheres  are  classified  as 
being  in  one  of  the  following  three  regions: 

Region  1:  (iZ,)  includes  only  the  Q  =  0  or  center  sphere.  It  combines  the  Z*  and 
A Zn  terms  which  have  the  largest  and  most  important  contribution  to  the 
impedance  sum,  and  its  contribution  is  evaluated  exactly. 

Region  2:  (A2)  includes  the  first  several  columns  of  spheres  a4jacent  to  the  x-axis. 
As  illustrated  in  Figure  H.3,  Aj  includes  spheres  in  the  rectangular  cylinder: 

-Nx<kM<N„  -N,  <  jv  <  N„  -oo  <  *,  <  oo,  (H.21) 

where  Ny  and  Nt  are  typically  on  the  order  of  8,  and  it  is  understood  that  the 
Q  =  0  or  R\  sphere  is  omitted.  For  R2 ,  two  approximations  are  made.  First, 
the  fields  of  are  approximated  by  the  fields  of  an  infinitesimal  dipole  of  the 
same  unit  current  moment.  Second,  the  weighting  function  W]  is  approximated 
by  a  Dirac  delta  function  of  the  same  current  moment  located  at  the  center  of 
the  center  sphere.  Basically,  the  Galerldn  weighting  function  is  being  replaced 
by  point  matching. 

Region  3:  (As)  includes  all  remaining  spheres  external  to  R2.  As  discussed  above, 
it  is  practically  impossible  to  directly  sum  the  impedances  associated  with  the 
As  spheres  to  convergence.  This  is  especially  true  if  he  is  complex,  since  in 
this  case  the  C**  in  Equation  (H.20)  will  exponentially  grow  as  (jy,  k,)  increase 
in  magnitude.  To  avoid  this  problem,  the  currents  in  the  A3  spheres  will  be 
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approximated  by  a  continuous  current  whose  field  can  be  easily  found  in  closed 
form. 

As  in  J?3,  the  fields  of  the  expansion  functions  in  A3  are  approximated  by  the 
fields  of  infinitesimal  dipoles  of  the  same  unit  current  moment,  and  the  Galerkin 
weighting  is  approximated  by  point  matching.  Since  the  separation  between  the 
spheres  is  electrically  very  small,  and  since  the  distance  from  a  R3  sphere  to  the 
origin  (field  point)  is  at  least  several  lattice  spadngs,  the  infinitesimal  current 
elements  in  A3  can  be  approximated  by  a  continuous  volume  current  density, 
J c,  with  the  same  current  moment  per  lattice  cell.  As  described  below,  the 
fields  of  J c,  which  exists  only  in  A3,  are  evaluated  as  the  fields  of  Jc  in  all 
space  minus  the  fields  of  Jc  in  Ai  +  Ri-  The  fields  of  Jc  in  all  space  are 
known  in  simple  closed  form,  and  the  fields  of  Jc  in  Aj  +  R2  are  evaluated  by  a 
relatively  fast  double  nunferical  integration  over  Ri  +  R2-  Essentially,  the  triple 
sum  to  infinity  is  replaced  by  this  finite  double  numerical  integration. 

Combining  the  contributions  from  the  three  regions  defined  above,  the  single  term 
impedance  matrix  element  can  be  written  as 

Zu  +  AZ11  «  Zri  +  Zr2  +  Zr 3  where 

Zr\  =  (#®  +  AZn) 

Na  oo  > 

z«  =  -  E  £  09  E  ft.(AQ) 

k,=-N,  jf=-N9  i«=— 00 

ZR3  =  -£cs,  (H.22) 

where  f?«*(AQ)  is  the  x  component  of  the  electric  field  of  the  infinitesimal  dipole 
located  at  AQ  approximating  J?,  and  Ec*  is  the  x  component  the  electric  field  of 
the  continuous  current  distribution  Jc  in  A3.  In  both  cases,  the  currents  axe  radiating 
in  the  homogeneous  host  medium  and  the  field  point  is  at  the  origin.  The  prime  on 
the  triple  summation  indicates  that  the  Q  =  0  center  term  is  omitted. 
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Figure  H.3:  Top  and  side  views  showing  the  different  current  regions  and  approxima 
tions. 
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The  Region  1  Term 


The  self  impedance  contribution,  Z/n,  is  the  contribution  from  Jf ,  the  expansion 
function  in  the  Q  =  0  center  sphere,  and  includes  both  Z®  and  A Zn.  Inside  the 
center  sphere,  E® ,  the  field  of  Jf ,  is  found  by  subtracting  the  incident  field  from 
the  total  field  in  the  eigenfunction  solution  [37].  The  total  field  is  approximated  by 
the  lowest  order  eigenfunction  term,  but  the  incident  field  has  an  infinite  number 
of  terms  of  slowly  decreasing  strength.  However,  due  to  the  orthogonality  of  the 
wavefunctions  when  integrated  over  the  center  sphere’s  volume  V®,  only  the  lowest 
order  term  contributes  in  the  integration  for  Z® .  Carrying  out  this  operation,  and 
also  performing  the  integrations  of  Equations  (H.9)  and  (H.10),  Zm  is  evaluated 
exactly  as 


AnI  — 


zm  =  (2*  +  A  Zn)  =  -  b\BM] 

St  ^  &i  sin  pi  cos  po  —  *o  cos  pi  sin  po  sinpisinpo 


3fcofci  t 


*?-*g 


*oPi 


Bint  — 


8ir 

3*o*i 


*osinpi  cospo  —  *i  cospi  sinpo  cos  px  cos  po 


*?-*$ 


*0Pl 


*osinpi  cospo  +  *1  co«Pi  nnpo  sinp1sinpo 


K \p\ 


k)PoPl 


(H.23) 


Region  2  Terms 


For  a  sphere  in  R2,  the  expansion  function  J?  is  approximated  as  an  infinitesimal 
dipole  of  the  same  unit  current  moment  located  at  the  center  of  sphere  Q.  Using  the 
well  known  expressions  for  the  electric  field  of  an  infinitesimal  dipole  [39],  we  obtain 

Zm  =  -  £  E  MAQ),  (H.24) 

kt=-N,  *m =-00 


where 


E„( AQ)  =  ^  Inn 


0  (’ + ;*or  +  ov)0 


where  a  is  the  angle  between  AQ  and  the  n-axis,  and  r  =  |AQj. 


(H.25) 
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In  evaluating  the  innermost  summation  of  Equation  (H.24),  the  field  Egx( AQ)  is 
summed  for  —  oo  <  *x  <  oo.  This  inner  sum  is  the  contribution  from  the  column  of 
dipoles  indexed  by  (jw,h«),  and  the  details  of  this  evaluation  are  given  below. 

Column  of  Dipoles  Summation 

This  section  evaluates  the  contribution  to  Zjn  from  a  column  of  infinitesimal  dipole 
current  elements  indexed  by  (jv,  kt),  i.e.,  it  presents  the  evaluation  of 

Z»k'  =  f;  ESx( AQ),  (H.26) 

i#— — oo 

where  Egx (AQ)  is  given  by  Equation  (H.25)  and  the '  indicates  that  the  ix  =  0  term  is 
omitted  if  jv  =  kM  =  0.  Experience  shows  that  Equation  (H.26)  is  a  slowly  converging 
sum,  and  simply  summing  to  convergence  significantly  increases  the  total  CPU  time. 

7 

The  method  presented  here  begins  by  directly  summing  terms  for  which  |tx|  is  not 
very  large.  Since  the  spading  dx  between  dipoles  is  electrically  small,  for  large  |*x|,  the 
discrete  array  of  infinitesimal  dipole  currents  are  approximated  as  a  continuous  line 
current  distribution  of  the  same  current  moment  per  unit  length.  The  field  of  this 
continuous  current  distribution  is  determined  by  an  asymptotic  integration  extending 
to  infinity,  which  can  be  approximated  in  closed  form. 

From  symmetry,  the  mutual  impedances  of  the  spheres  at  x  =  ±ixdx  contribute 
equally  to  Equation  (H.26),  and  thus  the  sum  in  Equation  (H.26)  can  be  reduced  to 
only  positive  values  of  ix.  Z’*k‘  can  now  be  written  as 

Z **•  =  ESx(yjvdy  +  ikxdM)  +  2  £  Egx( AQ)  +  ~  [°°  Etx( AQ)  dx ,  (H.27) 

<.=i  ®* Jx* 

where  Nx  is  the  number  of  terms  (dipoles)  that  are  directly  summed,  and  xa  = 
(JV.  + 1)  dx  is  the  x  coordinate  of  the  bottom  of  the  line  current  approximating  the 
dipoles  Nx  <  ix  <  oo.  In  Equation  (H.27),  the  first  term  is  the  sx  =  0  term.  It  is 
understood  that  this  term  is  omitted  for  the  jv  =  kt  =  0  column  of  dipoles  since  it  is 
the  R\  self  impedance  term  that  is  evaluated  separately. 

The  tail  end  contribution  to  Z**k‘,  given  as  the  integral  in  Equation  (H.27),  is 
evaluated  asymptotically.  Assuming  Nx  is  chosen  large  enough  that  xa  »  p  = 
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sjijydy)2  -f  (ksdx)2,  the  radial  distance  r  =  |AQ|  is  approximated  simply  as  z.  Thus, 
cos  a  «  1  and  sin  a  «  p/z.  Making  these  substitutions  into  Equation  (H.25),  the 
integral  contribution  becomes 


j r  £<,(AQ)<fe 

*  4*  [U)*4  o*o)>) i>+ 

(H.28) 

where  the  integrals  Jp  for  p  = 

1, 2, . . . ,  5  are  defined  as 

fOO 

J'=£  „ 

(H.29) 

The  p  =  1  integral  can  be  evaluated  by  replacing  the  exponential  by  its  Taylor  series 
expansion  and  integrating  term  by  term,  yielding 


h  = 


e-jko*. 

jk oz, 


«*  l 


1  - 


1  * 

+ 


6 


24 


jk oZ«  (jk oZa)a  (jkoXa)3  (JkoXa)* 

The  remaining  integrals  are  evaluated  through  the  recursion  relationship 

1 


(H.30) 


Ip  p-il*r 


g-jk 0*m 


(H.31) 


Region  3  Terms 


For  jR3,  the  total  currents  in  the  spheres  are  approximated  as  the  continuous  volume 
current  Jc-  J c  includes  the  e-,lt*'B’  variation  of  the  plane  wave  in  Equation  (H.l), 


i.e., 

3C  =  xJc  =  jtC2e-ik'{vco,'t’+*,in*)  in  Ra,  (H.32) 


where  C2  is  a  normalisation  constant  such  that  3c  has  unit  current  moment  in  each 
lattice  cell. 

If  Jc  existed  throughout  all  space,  instead  of  only  in  Jfe,  then  its  electric  field 
could  be  found  directly  from  Maxwell’s  equations  as  simply 


xC2 


jvp 0 
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To  find  the  field  of  3  c  in  Rs>  it  is  necessary  to  subtract  the  contribution  from  that 
portion  in  R2  and  R2,  yielding 

Zm  =  ~Ec,  =  -  -p-  f  JcB^(kop)dy  dx,  (H.33) 

*0  4W€q 

where  the  integration  is  over  the  region  R\  +  R2  in  the  yz  plane  (where  3c  is  zero), 
p  =  y/y 2  4-  z7  and  H^\k  •>)  is  the  cylindrical  Hankel  function  of  the  second  type. 

H.3  Numerical  Results 

This  section  presents  numerical  results  based  upon  the  above  PMM  analysis  of  an 
artificial  dielectric  composed  of  small  dielectric  spheres.  The  data  will  include  the 
root  of  ihe  fundamental  equation,  as  well  as  the  relative  effective  permittivity  and 
loss  tangent  for  different  a.  ^ficial  dielectric  compositions.  In  all  cases,  the  direction  of 
propagation  is  u*  =  z,  the  polarization  of  the  electric  field  is  x,  and  the  host  medium 
is  free  space. 

The  first  set  of  data  illustrates  the  root  of  the  fundamental  Equation  (H.12)  for 
the  effective  wavenumber  in  the  complex  K  plane.  The  geometry  of  the  artificial 
dielectric  consists  of  spheres  of  radius  0.0025Ao  spaced  in  a  cubic  array  where  dx  = 
dy  =  dx  =  O.OlAo-  The  spheres  have  relative  permittivity  eir  —  10  and  loss  tangent 
tan  =  1.  Figure  H.4  shows  the  magnitude  of  Zn(fce)  +  AZn  (i.e.,  the  magnitude 
of  the  determinant  of  the  order  N  =  1  matrix)  along  several  lines  in  the  complex 
k,.  plane  parallel  to  the  Re(fce)  axis,  as  illustrated  in  the  insert  to  the  figure.  The 
root  is  the  value  of  K  such  that  |ZU  +  AZn|  =  0.  From  Figure  H.4  this  occurs  at 
Ae/fco  ~  1.086  — *0.0125.  Using  Equation  (H.2),  this  corresponds  to  a  relative  effective 
permittivity  of  s«r  =  1.18  with  effective  loss  tangent  tan  Se  =  0.023.  Note  that  the 
effective  permittivity  of  the  artificial  dielectric  is  between  that  of  the  host  medium 
and  the  dielectric  spheres. 

The  insert  in  Figure  H.5  shows  an  artificial  dielectric  in  which  the  dielectric  spheres 
have  relative  permittivity  eir  =  4  with  loss  tangent  tan  6i  =  1-  The  relative  effective 
permittivity  and  effective  loss  tangent  are  plotted  for  sphere  radii  of  a  =  0.002,  0.01 
and  0.02Ao,  as  a  function  of  cubic  array  lattice  spadngs  <f,  =  dy  =  dM  varying  km  2a 
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Figure  H.5:  Effective  permittivity  versus  array  lattice  spacing,  for  a  3D  array  of 
dielectric  spheres. 


to  O.lAo-  As  the  lattice  spacing  increases,  the  fraction  of  the  volume  occupied  by  the 
spheres  decreases,  and  the  complex  permittivity  of  the  artificial  dielectric  approaches 
that  of  the  host  medium. 

The  insert  in  Figure  H.6  shows  an  artificial  dielectric  in  which  the  dielectric 
spheres  have  radius  a  =  0.0025Ao,  and  are  arranged  in  a  cubic  lattice  of  spacing 
dx  =  dy  =  dx  =  O.OlAo-  The  relative  effective  permittivity  and  loss  tangent  axe  shown 
as  a  function  of  sphere  loss  tangent,  tan  Si,  varying  from  0  to  10,  for  sphere  relative 
permittivities  of  €ir  =  3,  10  and  30.  Note  that  the  relative  effective  permittivity 
increases,  whereas  the  effective  loss  tangent  decreases,  with  increasing  eir.  When 
tan  Si  =  0,  the  spheres  are  lossless,  and  thus  the  artificial  dielectric  is  also  lossless 
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Figure  H.6:  Effective  permittivity  versus  sphere  loss  tangent,  for  &  3D  array  of  di¬ 
electric  spheres. 

with  tan  6e  =  0.  As  tan  6i  increases,  tan  6e  increases,  reaches  a  peak,  and  then  falls 
to  zero  as  tan  6e  —*  oo. 

H.4  Summary 

This  appendix  has  presented  an  integral  equation  and  PMM  solution  to  determine 
the  effective  permittivity  of  an  artificial  dielectric  composed  of  a  3D  periodic  array 
of  homogeneous  dielectric  spheres  in  a  homogeneous  host  medium.  The  effective  per¬ 
mittivity  is  determined  by  finding  the  complex  wavenumber,  K,  for  a  plane  wave 
propagating  in  the  artificial  medium.  The  method  can  compute  the  effective  permit- 
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tivity  of  the  artificial  dielectric  u  a  function  of  frequency,  direction  of  propagation, 
and  the  size,  material  composition,  and  density  of  the  dielectric  spheres. 
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